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ABSTRACT 



A preprocessor is designed to extract a set of features 
that enhance natural clustering by removing extraneous 
information. The design removes time shift and scale 
dependence by taking advantage of invariant properties of a 
Fourier transform followed by a Mellin transform. The 
preprocessor is realized using an FFT and a tfellin 
transform with a conventional error correction term. The 
error term proves to be indeterminate, but the error's 
bound is identified as the envelope for Beilin correction 
terms. Properties of the Beilin transform are employed to 
modify the signal so that the error correcting is no longer 
required. The resulting algorithms are tested with 
variously scaled inputs for which closed form solutions are 
known. With a verified modification in place, the 
preprocessor produces features that are invariant to 
shifting and scaling, while retaining enough information to 
classify canonic shapes. A method of improving performance 
is introduced. 
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I . INTRODUCTION AND BACKGROUND 



Pattern classification is the assignment of a physical 
object or event to one of several prespecified catagories 
and is the result of an incomplete theory of perception. 
Although many transducers are available for converting 
light, sound, temperature, reflected radar signals, etc., 
to electrical signals, the ability of machines to perceive 
cr to recognize their environment remains very limited. In 
the structured world of communcation s engineering, signals 
are designed to be detectable and differentiable. A much 
more difficult problem presents itself when sensing an 
environment through a transducer and recognizing or even 
classifying the elements of that environment on the sensed 
characteristics of the transducer's electrical output. 
Pattern recognition can be considered a complex 
communications problem (for example, attempting to teach a 
machine to decode signals encoded by nature). It is 



possible to alter the 


t ra nsducer ' s 


output 


to facilitate 


object classification, but 


determining how to 


alter 


that 


output is not a simple 


task. The 


mai n 


elememt s 


o f a 


classification system is 


shown in 


figure 


1 [1]. 


The 



transducer senses, actively or passively, a set of 
characteristics belonging to the object. These 
characteristics can never be a complete description of an 
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TRANSDUCER 




A Classification System 
Figure 1 
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object, but represent, hopefully, enough information to 
classify the object as belonging to one of a number of 
classes. For instance, temperature is a characteristic of 
the object and a class, but this feature is of little value 
unless it differs in some way from objects belonging to 
other classes, while the set must include characteristics 
that are common among that class. The preprocessor (or 
feature extractor) aims to reduce the data by measuring or 
quantifying certain properties that distinguish the sensed 
object as belonging to one class and not to others. This 
can be done by discerning key features or using the imput 
to generate another set of features optimized by some rule. 
The values of each of these features is then passed to the 
classifier, which evaluates these features to assign the 
object to a class. 

With varied success, machine pattern classification has 
been applied to a large range of problems/disciplines. 
Fields where it is particularly common include optical 
imagery, acoustic signal processing, radiology, radio 
astronomy, and electronic warfare, to name a few. Work in 
many of these fields was reviewed during the development of 
this thesis and the results derived and demonstrated here 
may, in turn, be applicable to the field in general. This 
effort has been directed toward designing a preprocessor to 
produce a set of features that are invariant to information 
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known to be superfluous to classification, but that retain 
enough information to classify an object. The object has 
been sensed by a transducer and has been represented as an 
empirically derived, univariate time series. Such a series 
would be the form of data available from range only radar 
return which is specifically what the preprocessor is 
designed to handle. Returning to figure 1, the object has 
an infinite set of characteristics (here portrayed as an 
infinite series of discrete values). The 
transducer/receiver has collected some characteristics of 
the target objective in the presence of noise. This 
information is relayed as a set of discrete samples (hi) 
from a band limited signal. The preprocessor is designed to 
determine and code revelent information (Hj) for the 
classifier. If this task was done well, classification 
becomes a trivial problem. On the other hand, if the 
classifier becomes ideal (capable of resolving an infinite 
number of characteristics in noise) the preprocessor design 
begins to look like a wire. The distinction between the 
preprocessor and the classifier is arbitrary from an 
analytical point of view. When designing a classification 
system functionally, a difference is enforced. The 
classifier has little concern for how the features are 
developed, but seeks to efficiently use those provided to 
guess the class of the target object. The preprocessor is 
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problem dependent, needing to produce an optiiral set of 
features, Ej , from the sensed data hi. 

A. FEATURE EXTRACTION 



For the purposes of this paper, two generic approaches 
to feature extraction are defined. The first, a 
classification approach, was described above. The second, a 
descriptive approach, tries to define the object in terms 
of the objects' structural features. This system might 
recognize a car, for example, by breaking up the visual 
picture into canonic shapes, and comparing this to 
previously specified canonic class models. The perceived 
structure of the physical object is maintained and should 
reflect the structure of the object itself. This approach 
could be robust to temporary changes in the object itself. 
In the car example, knowing that at one end of the car the 
trunk can be opened or closed allows the device to take 
this factor into account. Another important advantage to 
descriptive techniques is that the class characteristics 
may be entered or specified without collecting actual 
transducer generated data to train the machine. 
Unfortunately, the problem of designing a machine to 
analyse a visual scene to produce a structural description 
has proved to be quite difficult. Object description from a 
univariate time series is even more difficult, and if the 
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radiation sensed by the transducer is not frorr the visual 
spectrum, the task rapidly approaches the impossible. For 

i 

these reasons the approach taken was the classification 
approach (to reduce the signal to a set of orthogonal 
features that do not uniquely reflect the structure of the 
object, but do retain sufficient information to classify 
the ob jec t ) . 

This paper excludes a detailed description of the 
transducer specification. The problem of the classifier 
itself is viewed as one of partitioning the feature space 
(Hj) into regions? one region for each category. Ideally, 
this partitioning should be arranged so that none of the 
decisions are ever wrong. When this cannot be realized, at 
least the probability of error should be minimized or the 
average cost of errors minimized. The problem is one within 
statistical decision theory. Knowledge of the object 
classes (the transducer and the classifier) are required to 
design the preprocessor, which is the topic of this thesis. 
The preprocessor designed and built here generates features 
from a range only radar video signal. These features are 
used by a general Bayesian learning classifier. The 
supervised learning general Bayesian classifier approaches 
the problem by taking a series of incoming sets of features 
labeled as to their class. From the data, an a posteriori 
density is computed. Sach successive set of training data 
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is used to refine the densities' statistics. Vhen the 
classifier has been trained on N classes, the features are 
rrodified to separate the class volumes in an optical way 
and to reduce the nurrber of features to one less than the 
number of possible classes. The feature vectors of class 
members are clustered about a simplex point and likely 
boundaries are set up allowing classification of the object 
as belonging to one of the classes, or of an unkown class. 
An N simplex is a collection of N points in ( N— 1 ) space 
where the distance between any two of the points is equal 
to the distance between any other two. Thus a three class 
problem transformed into a three simplex in a two 
dimensional plane produces clustering of the three class's 
about the vertices of an equilateral triangle. The simplex 
coordinates are the reduced feature vectors, generalized 
from the training data. In a controlled, simple problem the 
classifier works well, but when encountering real problems 
one class's feature space will intertwine another's, making 
it much more difficult to obtain separation in a meaningful 
way. The goal of the preprocessor is to present key 
features that determine class for subsequent optimization 
by the classifier. 
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B. FOUR IIR - MELLIN PREPROCESSING 



Common to all univariate, time series classification 
problems are several variables that interfere with the 
recognition process. Assuming discrete data processing is 
used, these are addressed in the following order; 
windowing, framing, scaling, sampling rate, quantization 
noise, and sufficient information. For real processing, the 
imput waveform is not sampled for all time. It is sampled 
for a period of time. This windowing of the data corrupts 
the resulting spectrum in two ways [2]. First it introduces 
a periodicity (the inverse of the window length) and 
resulting aliasing to the otherwise infinite spectrum, and 
further distorts the spectrum by a convolution with the 
spectrum of the window itself. Both of these effects will 
color all of the data in the same way and so can be 
accounted for by deterministic methods.' Framing can be 
considered characteristic of poor synchronization, whereby 
the pattern of concern is not position stable with repect 
to the window as shown in figure 2a. It seems that even in 
human optical recognition, the eye tends to center the 
pattern prior to recognition, unless trained otherwise. 
Scaling is that property whereby the object field may vary 
in scale or aspect angle, in one dimension or several 
dimensions as in figure 2b. Before the analog signal is 
sampled, it must be fed through a sharp low pass filter. 
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Shifting and Scaling Variation 
Figure 2 
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because no higher frequency noise can be present without 
being folded onto the valid data. Quantization noise, due 
to the requirement to round off each sample to some 
discrete level, is treated the same as round off error in 
numerical processing [3]. In all of these problems 
discussed to this point, the effect of this processing is 
to mask the actual feature vectors, making the 
classification system less sensitive to valid pattern 
characteristics. In all recognition problems, it is assumed 
that there is sufficient information present for a pattern 
to be detected and identified by the system. This means 
that there is sufficient variability between the classes, 
but sufficient similarity between those patterns of a 
single class to classify each pattern in terms of those 
classes . 



A verifiable goal of this preprocessor is to produce a 
set of features that are invariant to shifting and scaling 
changes. An approach, figure 3, has been proposed and used 
[4-10] . The preprocessor consists of putting the sampled 
data successively through two transforms, a fast Fourier 
transform and a discrete Mellin transform. Integral 
transforms, the Fourier and Mellin included, develop 
naturally from the solution of simple problems in potential 
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theory [11-13]. The Fourier integral transform of the 
waveform h(t), 

H(f) = ^ ‘id) k(f, t)4t (1) 

~fc> 

where 5=exp [-2 iff t] , is a principle analytical tool in such 
diverse fields as linear systems, optics, probability 
theory, quantum physics, and signal analysis [13]. Its 
purpose in the preprocessor is twofold, but relies on a 
single characteristic. The magnitude of the Fourier 
transform is invariant to shifting, h(t-a). 

I I = / H( s)l ( 2 ) 

This characteristic removes the effects of framing 
inaccuracies, and also permits the averaging of successive 
looks or pulses of data to improve feature resolution, but 
removes much of the information about a signal's structure 
as discussed later. A discrete Fourier transform, 

( t ) -> 4 ) tyy\ - 6 , l j "Zy • . « i M — I 

H ( f) —> H </*> St) * o , 1 , 2 j j N ~ \ ( 3 ) 

M-l 

n V / . • 

H(rr>) « / ■ € ( 4 ) 

~0 

has an equivalent identity, but it is only exact for shifts 
of integer values. 
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(5) 



I J* C A<*-+) ] I 

- I Hc~> -|H(-)J 

Shifts of other than integer values result in errors that 
depend not only on the shift, hut on qualities of the 
sampled waveform itself, h(t). 

The Mellin transform is an integral transform with the 
kernel, K. 

K - t ( 6 ) 

h ui « j'ia; t'-'wt (7) 

is the Mellin transform with respect to the complex 
parameter s=r-j2-f>m. Several simple substitutions relate 
this to more common analytical tools. Exponentially warping 
t=exp[x] changes the appearance of the integral to what is 
often called a double sided Laplace transform [14], 

Hcsj = C £(&*) ^ 

J -eb 

The transform is invariant to t domain scaling when taken 
with respect to the imaginary part alone, 

= f "°X '(*) & (9) 

-ft? 
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Equation (S) is recognized as the Fourier intergral of an 
exponentially distorted waveforir h'(x). The modulus of the 
expression on the right is the magnitude of the Fourier 
transform of the exponentially distorted time function. The 
property to be exploited in a Fourier-Mellin (FM) 
preprocessor is that the modulus of the transform in s, is 
invariant to t-scaling. Given a time waveform h(t), its 
Beilin transform H(s) is given as equation (7). Scaling the 
entire t-domain by k, 



■mtA <*/*)]* (10) 

Letting r=t/k f and remembering that s is imaginery, 

•k J A(.r) 

jjk* MCs) I a I HCS) I (11) 



Much effor 
digitally 
is devoted 
di screte 
properties 
treatment 
concepts 
content. 



t and detail is spent implementing equation (7) 
in Chapters II and III. The rest of this chapter 
to providing some required background on the 
versions of the Fourier transform and some 
upon which the FM preprocessor depends. The 
here is brief, being mainly a review of basic FFT 
and as such may be skipped without loss of 



A discrete 
efficient and 



Fourier transform is not computationally 
so leads to impractically long processing 
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times. The fast Fourier transform ( FFT ) efficiently 
computes the discrete transform and is used in the 
preprocessors "built for this thesis. Other properities of 
the FFT should he presented before getting into the 
detailed design of the preprocessor. The first and last are 
concerned with symmetry. If h(m) is real, as in the case of 
the sampled data, then the frequency spectrum of that data 
is even, 'Ee(n) 1 ='Ee(-n) ! . B(n) has a real part and an 
imaginary part, F.e(n) and Im(n), while h(m) has an even 
he (rr )=he (-m ) , and an odd ho (m) =-ho (-m) part. 



The odd part of h(n) times the cosine kernel summation, and 
the even part of h(n) times the sine kernel summation are 
both zero. E(n) can now be seen to have an even and odd 
part. Taking the magnitude of an odd function makes it 
even, proving that the Fourier transform of a real series 
is a spectrum whose magnitude is symmetric about f=0. This 
is of course true for both the integral and discrete 
Fourier transforms. The second line of symmetry is an 
effect of discretizing the signal and its spectrum for a 
discrete transform as evinced by its theoretical 
development. Before proceeding though, the convolution 
theorem is required. 





( 12 ) 
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The convolution of the two functions h(t) and g(t) is 
defined as the farriliar integral, 



The relationship between convolution and the Fourier 
integral is very important to modern analysis and 
contributes to making the Fourier transform a key analytic 
tool. The convolution theorem states that if h(t) has a 
Fourier transform F(f) and g(t) has the Fourier transform 
G(f), then h(t)#g(t) has the Fourier transform H(f)G(f). 



Proving this, the Fourier integral is used directly on the 
convolution integral. 




( 13 ) 





( 14 ) 





( 15 ) 



From the shifting theorem already presented, 



Y(f>* C -tr 





( 16 ) 
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And finally, 



Y(-f) ' J [ j c-t) *J!c*)j= Cjct)H(-f) 



(i?) 



It can be shown similarly that. 



¥ hco] = ^ (t)jict) 



( 18 ) 



Clearly, convolution in one domain is simple multiplication 
in the other domain. Although not needed for the pending 
development of the discrete Fourier transform, another 
important relationship, known as the correlation theorem, 
can be appropriately dealt with here. The correlation 
integral , 



has an operation with which it forms a Fourier pair as did 
the convolution-multiplication operations. This theorem can 
be established as before, 




(IS) 





( 20 ) 
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The term in brackets is the complex conjugate of H(f) and 
is denoted by E*(f) in the final form of the theorem below. 



3 [ J J? (f) 4 (r + t) oL'T'J - H*(t) 

"0O * 



( 21 ) 



We will now continue on to the discrete Fourier 
transform starting with a continuous waveform h(t). The 
waveform is sampled or multiplied by a string of delta 
func tions, s( t) . 



where capital delta is the sampling interval. The infinite 
sum is not realized and must be windowed, in this example, 
by w(t)=l for 0^t £(M-1)4 and zero elsewhere. So that 
now , 



Recalling the convolution theorem, the multiplications in 
the time domain correspond to convolutions in the frequency 
domain with the following results. H(f) is convolved with 
the window functions' spectrum and will have the apparent 
effect of introducing ripples because of the window's 
significant sidelobes. The rippling may be minimized by 
choosing a window function with small sidelobes at the cost 
of other, perhaps more acceptable, spectral degradation. 




( 22 ) 




( 23 ) 
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E(f) is convolved with the sampling function and 
S ( f ) =1 ( f-n/4, ) has made the spectrum periodic with respect 
to the interval 7=1/4. . The spectrum coming from a real 
waveform is first symmetric about f=0 as shown before. Now 
because of its periodocity the spectrum is symmetric also 
about f=F/2 (the Nyquist or folding frequency). One final 
step remains. The Fourier spectrum is also taken at 
discrete points 1/T apart. The result in the time domain is 
the convolution of the sampled, windowed signal with 
I(t-nT), which is a periodic signal with T as its period. 

The FFT algorithm used in the preprocessors developed in 
this thesis uses a Cooly-Tukey, base two algorithm [15]. 
This is documented where it is used in programs included in 
the Appendices. The algorithm uses N samples where N is two 
to an integer power. In the transformed domain, due to the 
symmetry shown above, there are N coefficients (only N/2 of 
which are unique as shown in figure 4). The original 
waveform must be band limited prior to sampling to minimize 
aliasing. The resulting frequency spectrum should approach 
zero at the folding frequency where up to half of the power 
can be aliasing noise. 

In the following two chapters, several different 
preprocessing algorithms are developed and their 
performance compared with canonic inputs. In Chapter IV the 



2 e 



problem of ship identification with range only radar is 
discussed briefly and one preprocessor is applied to 



several ship profiles. In the fifth and final 
conclusions are drawn from the work done here and 
efforts are recommended. 



chapter , 
follow-on 
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An FFT Spectrum 
Figure 4 
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II . DIGITAL ^SLLIN TRANSPORT 



3T EXPONENTIAL WAPPING 



In the first chapter, the modulus of the tfellin 
transform was shown to be invariant to scaling. A detailed, 
examination of the mechanics involved suggests an 
implementation that is widely used. The Mellin transform of 
a t-domain function h(t) is given in equations (?) and 
O'). 



Delta has been added, corresponding to the sample interval 
in the t-domain. Again it is noted that O') is a Fourier 
transform, where s=-j2mr. Solving the integral for the 
effect of a t-scaling by the factor k, 



Clearly, in the Fourier integral, the scaling factor k has 
become a shift for which the modulus of the transform is 
invariant. The exponential warp alone has transformed the 
scale factor into a shift. A prerequisite is that the 
t-domain signal has no shift. If there were a shift, it 
does not transform to a simple factor or shift in the 



His) - £ (*) 



(?) 




(s') 





(24) 



2S 



Mellin dorrain and so cannot subsequently be rerroved by 
taking the magnitude of a Fourier transform. 

Implementation of a discrete Mellin transform is as 
difficult as it is with the Laplace transform [13]. Once 
transformed, characteristics in the new spectrum are 
difficult to relate to the original signal characteristics. 
One hypothesis relating the two domains is generated by 
comparison to the Fourier transform. The power spectrum 
associated with the Fourier transform can be used to detect 
periodicities in the physical function, since the wave 
numbers at which sharp peaks of the spectrum occur give the 
wavelength of such periodicities. 3y analogy, the positions 
of the peaks in the spectrum associated with the Mellin 
transform are said to give the magnification or compression 
which will produce features in the physical domain. 
Further, this stretching and compressing is identified as 
periodic in nature [21]. This seems unlikely because the 
Mellin is invariant to magnif ication/compress ion and does 
not behave well (a scale factor that is a function of t and 
k(t) as seen in (24). The Fourier spectrum models the 
original signal by a set of weighted sinusoids of varying 
frequencies, and therefore, naturally display periodicity 
and is invariant to shift. The Mellin is also a weighted 
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suit of sinusoids, but whose magnitudes are inversely 



weighted by t. 



His) - 




Aw 



_ 3 “t 

e * 






(25) 
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Values for h(t) for 0 t 1 are far more important to the 
sum than those beyond that point. This characteristic 
contributes to the difficulty of realizing discrete Laplace 
and Mellin transforms and is a major topic covered in this 
chapter . 

The numerical approximation of the Fourier-Mel 1 in (Fi^) 
transformation by exponential warping is functionally 
diagrammed in figure 5. A FORTRAN program using the 
algorithm described in this chapter is included in Appendix 
A. Referring to figure 5, the input samples are from a 
pulse whose duration is finite and less than that of the 
sample window. For this case, no spectral distortion is 
experienced by filling zeros behind the sampled data. The 
only effect of the zero filling is to add spectral 
resolution in the frequency domain. Once through the first 
FFT block, and the magnitude taken, Ef is symmetric about 
zero and at the folding frequency. Thus, for N time samples 
and filled zeros, there will be N/2 unique spectral 
coefficients. This unique spectrum is interpolated, 

resampled as a warped signal and fed to the final FFT 
block. A correction is added, and the modulus is taken. The 
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F M FEATURES 



Fourier-Mellin by Exponential Warping 
Figure 5 
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resulting FM features are invariant to shifting and scaling 
in the time domain. The FIT block was covered in sufficient 
detail in the preceding chapter. Some effort will be spent 
in discussing the warping itself and the need for, and the 
development of, a zero point correction. 

A. ALGORITHM DEVELOPMENT 

This section uses its own notation to address the 
requisite exponential sampling. The series to be 
transformed is h( f ) . Its Mellin transform is E(m) where m 
is the Mellin frequency in s=-j2^m. The transform, is, 

fk'v-O = f £(■£)£ (2 6) 

Letting f=Fexp[x] as before, where F is added corresponding 
to the sample interval 

H(^) - J ACfe.*-) (27) 

The need is to evaluate M equally spaced samples in x, at 

O t X, XX, . . . , < M- I) x (28 ) 

while the data consists of N equally spaced samples in f, 

0 , F, T-F, • < ( N-0 F ( 2S ) 

Assuming that F is a small enough interval to properly 
characterize h(f), it is sometimes advisable [9,10] to 
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choose the exponential sampling interval (X) such that the 
largest intersample spacing in h(Fexp[x]) is equal to kF 
where k=l . The other set of conditions used to uniquely 
specify the new samples are that f=F and x=0 will he the 
lowest limit for interpolation, while (N-l)F and (P-l)X are 
equated as the upper limit. The first requirement 
constrains the choice of M by 

A = e (M ' z)X ( 30 ) 

Meeting the second condition, the end points in each 
sampled series are equated yielding, 

(n-i)* e CM ' lJX (3i) 



Substituting (31) into (30) while applying the exponential 
series approximation gives. 






(32) 



N - / 



One more substitution, (32) back into (31) sets up the 
desired result where M is now specified to exponentially 
sample from f=F to f=(N-l)F with kF being the largest 
interval between samples. 



M 



A. c aw-/) 

* 



(33) 
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As N gets larger, M explodes. If N=l€ then P=41 , if N=32 
then M=106, if N=64 then M=261, and so on. This Strict 



requirement can be compromised depending upon the 
application. The other extreme [18] is the criterion that 
requires the smallest interval between samples to equal the 
interval between uniformly spaced Nyquist sampling, with 
the intervals increasing exponentially thereafter. The 
specification permits analyzing frequencies approaching the 
largest which can be analyzed with uniform spacing. This 
greatly reduces the number of samples and decreases the 
required processing time, but is limited in application. 
Two factors mitigate the stringent requirements imposed by 
(33) where k=l . First, using N unique, uniform samples 
yields N/2 unique, uniform samples in the FF? domain. A 
related consideration is that the values of the spectral 
components approach zero at the folding frequency because 
the original signal has been band limited and some over 
sampling is normally recommended. Secondly, the inverse f 
weighting apparent in equation (25) 



attaches a decreased importance to h(f) near the folding 
frequency fn. These two effects combine to permit a much 
lower sampling rate once the modulus of the FFT has been 
taken. In spite of this, the stiffer rule (33) is used for 
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this work to generate the test FM dorrain possible. Whatever 
rule is used, once the exponential sample points have oeen 
computed for an FM preprocessor, they needn't be 
recomputed, but can be stored for rapid access during the 
interpolation. Some interpolation must be performed to 
approximate the spectral values of the new sample points. 
Third or even forth order Lagrange polynomials have been 
recommended and used for this purpose with apparent 
success [9,10,19]. The advantage of the Lagrange technique 
is that its notation is particularly compact, consisting of 
simple summations and repeated products that lend 
themselves to digital realization. Unfortunately, for the 
data sets used in this thesis, the third order algorithm, 
was observed to behave poorly, adding a ripple in regions 
of rapidly changing slope. That this might be the case was 
suggested by the advice that Lagrange is very good near the 
central data point when the order of the polynomial is 
known to be the same as the order of the approximation, 
otherwise it is test left alone [20,21]. In its place, a 
second order polynomial was used to interpolate the warped 
samples with results that were nearly indistinguishable 
from the actual waveform, as seen later in figure 8. 
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3. ZEPO POINT CORRECTION 



Another problem becomes apparent when the Mellin 
transform of h(f) is recalled, 

Hts) » f ; 'Af - (3 ^ (34 ) 

where s=-j27>m. The exponentially sampled waveform 
described above is applied to an FFT block. As f approaches 
the folding frequency, h(f) tends to zero. Unf or tuna tely , 
as f approaches zero, the value of h(f) is not zero. In 
fact it is frequently rather high. To make matters even 
more interesting, the left integral in equation (34^ 
clearly shows that the closer to zero f gets, the more 
important h(f) becomes to the integral. Several solutions 
to this problem are considered below. 

One practical, simple approach is to set the DC (ie, 
f=0) term of the FFT to zero. The effect is nothing more 
than removing a DC level back in the signal domain, but 
Mellin transformation into the FM domain leaves the 
spectrum dependent upon the scale factor k [22]. Setting 
the f=0 coefficient to zero corresponds to setting h(f) to 
zero for 0 f<l, where the unity upper limit is chosen 



3 ? 



without loss of generality. The resulting transform of a 
scaled signal domain h(f/k) is 






(35) 



which is obviously dependent upon k. The closer k is to 
unity, the smaller the effect. By increasing the f spectral 
resolution, the error can be reduced. The error may be 
insignificant for many applications [5-8], but the 
technique should be used with care. 

The first solution has highlighted the need for a zero 
point correction. Another common solution [10-11] is 
developed by breaking the integral up as before. Again 
using unity as the upper limit of the left integral while 
remaining general in application, 



Two assumptions are made to get the correction term. First, 
that h(f) remains a constant h(0) over the interval of f 
between zero and one. Second, and not as easily accepted, 



Equation (37) pretty well shows why this assumption is 
suspect, but playing along for the momment, the question is 







(36) 




(37) 
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detailed look, 
factor becomes, 







Accepting the 



(38) 



Because H(s) is a complex function, Z(s) must be applied 
(added to the imaginary part of the succeeding Fourier 
transform) before the magnitude is take to remove scaling 
dependence. This correction is specifically derived for use 
with a continuous Fourier transform such as the optical 
Fourier in imaging systems with the added stipulation that 
h(f) be nearly constant over the range 0<f<k where k is the 
largest scale factor expected. If an FFT is employed to 
make this final transform, another correction should be 
applied as shown by Zwicke and Kiss [11] below. This 
correction factor differs from the first due to an 
invariant property of the FFT. The FFT of two unit step 
functions that vary only in scale are balanced. 



IZj ♦ t | 



(3S) 



where m and p are arbitrary integers greater the zero and 
less that I*. Successive FFT coefficients are summed, and 
the average value of the contributing terms taken resulting 
in, 



c = 



— - J_ co-6 ( lr J: /m) 

z Z 



(40) 
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This is then multiplied by h(0) to arrive at the FFT Beilin 
correction factor. 



Z ppt C4)= ASlL (i- 



(41) 



When k/M is small the imaginary term dominates and the 
correction factor approaches that used in the continuous 
case (38). Most of the work done for the thesis on the 
exponential algorithm was done using the inappropriate zero 
point correction (38). Since its discovery was coincident 
with that of more powerful methods discussed in Chapter 
III, little data was taken using (41). 

To bound the error involed, an acceptable h(t) is 
defined, windowed and transformed using the Mellin 
integral. The window limit is then allowed to grow 
unbounded and the resulting expressions are interpreted as 
the error. It's been assumed that for 0<t<l, h(t)=h(0). 
Since the error resulting from the assumption in equation 
(37) arises from this interval alone, t>l is ignored for 
the time being. Warping the signal as before h(x)=h(0) for 
all x<0 and h(x)=0 for all t>0. Evaluating the integral to 
a finite window width (T) , 





( 42 ) 
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The term in brackets is the magnitude of the contribution 
from the region of integration, -T<x<0. Note that it 
involves a sin()/() term. The effect is to add a peak of 
( T ) h ( 0 ) at the origin. The size of the peak depends 
directly on the window border T. Letting T approach 
infinity raises the spike at <^>=0, with lesser peaks at 
^ =(2i+l ) tr/T, where i is an integer. Fach of the subpeaks 
has a magnitude of 2h ( 0) T/ ( ( 2i+l ) / tT* ) . Substituting uO into 
the second relation yields the envelope (in brackets) and 
phase as T tends to infinity. 



The error bound in brackets, does not depend on the sample 
rate, or the size of the window. Any approxima tion 
approaching zero will have the same bound. Although the 
magnitude of the error is in a convenient form, the phase 
is indeterminate. For the correction to be applied, the 
complex addition must occur prior to the modulus being 
taken. This cannot be done, leaving the error uncorrected, 
but is bounded by 2h(0)/*o for continuous and aoeriodic 
discrete Fourier transforms. For the FFT, equation (43) 
does not bound the error. The sum of 1/i does not converge 
as i tends to infinity. At any point on the FFT this sum is 
present due to the apparent folding. The error itself is 
not unbounded because phase differences in the sum of the 





(43) 
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errors at any point tray result In the envelopes adding 
destructively, reducing the actual error. Adding the 
envelopes is an unrealistic, worse case approach. The FFT 
correction (41), can not he compared to (43) for these 
reasons. The other error correction, using the assumption 
in equation (37), can he compared. The first, setting h(0) 
to zero, or just plain ignoring the 3<t<l interval have the 
error function hound in equation (43). Although only 
differring hy a factor of two in magnitude, the constant 
phase is arbitrary, and equivalent to setting T=0; that is, 
assuming h(0)=0 over 0<t<l. This was the very problem the 
correction was developed to remedy, hut is without effect. 
Equation (41), the correction for the FFT is not completely 
accepted hy this author, albeit no real empirical evidence 
have served to verify or dispute the claim. Suspicions are 
raised on two aspects. The indeterminate phase of the error 
(43) in the continuous case arises naturally from 
approaching the t=0 limit with the Mellin integral. This 
quality is conspicuous hy its absence in the FFT error 
correction. The FFT error correction is computed by summing 
the FFT coefficients in the complex plane. The average 
position of the resulting polynomial is the- error 
correction term. Fowever, the FFT coefficients of a finite 
duration signal are the values of the signal, evaluated at 
M evenly spaced points about the unit circle [3,23]. If 
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the sequence is a constant, the average value is the origin 
of the complex plane. 

The zero point corrections for the Mellin transform are 
unbounded at «^>=0. More time could have been spent 
determining the best applied correction for the specific 
case at hand, but direct methods are developed in the next 
chapter that obviate the need to employ the correction at 
all . 



C. TESTS AND RESULTS 



It is worth admitting at this point that the results 
using the exponentially warped algorithm to achieve a 
discrete Mellin transform have not been good. More recently 
developed techniques in Chapter III greatly surpass the 
results reported in this section. Although much of the 
theory used to improve performance in the following 
chapters could have been used here, this was a preliminary 
attempt that was later abandoned. The FM processor 
described in the previous section was built using FORTRAN. 
Appendix A is the documented program. This section will 
review the processing with actual plots of the signatures 
at different stages, discuss the required tests, introduce 
the testing approach, and finally intepret the results. The 
functional block diagram of the preprocessor, figure 5, is 
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near the beginning of this chapter and could be profitably 
reviewed. Figures 6 through 9 represent the step by step 
view of the signal processing, where the signal, a test 
shape, is shown in figure 6. The waveform, appears as an 
envelope, and is drawn with vertical lines indicating the 
sampled series. Figure 7 is a picture cf the FFT, in its 
sampled version showing its symmetries. Figure 8 is a very 
close approximation to the continuous periodic transform 
achieved by filling zeros to obtain the requisite spectral 
resolution. A "+" on the transform plot indicates an 
exponential sample point interpolated from the sixteen 
unique points in figure 7. The warped samples are sent 
through the FFT block once more with the result shown as 
figure 9. Heavy spectral coloring by the h(0)/ui correction 
factor is evident. Only the first half of the spectrum, 
from zero to the folding frequency, is valid. 

Complete scale invariance was never realized although 
its effect was greatly reduced. All the testing done on the 
exponential algorithm was an attempt at achieving and 
verifying shift and scale invariance. Much effort was spent 
structuring the tests to avoid the effects of processing 
noise so the actual algorithmic characteristics could be 
determined. Along the way, requirements arose to select a 
suitable interpolation polynomial, an optimal zero point 
correction and other system improvements. In all cases, the 
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Test Shapes 
Figure 6 
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Test FFT 
Figure 7 
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Test Exponential Sampling 
Figure 8 
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test procedure was to first select canonic test shapes. 
Squares and triangles were most frequently used, being 
shifted, scaled and combined to determine preprocessor 
characteristics. Scaling was most frequently by a factor of 
two or less. This corresponds to an aspect angle change of 
60 degrees from the unsealed case. For many tests, care was 
specifically taken to exactly scale a sampled series 
instead of the waveform envelope. The variation in input 
signal when this isn't done is evident when considering 
figure 10. For the envelope shown, the sampled series 
cannot reconstruct the same signal, and may result in 
feature space variation. Two feature qualities were 
monitored in each test to determine preprocessor 
performance? insensitivity to shifting and scaling, and the 
ability to differentiate between canonic classes. These 
qualities were measured by visual comparison of the 71* 
features, by computing the correlation coef f iecients and 
the mean squared error between the feature sets of 
differently scaled similar test shapes, and by computing 
the distribution of the error over the feature space. The 
latter test was an attempt to locate feature regions of 
class commonality, and regions of distinction between 
canonic classes. This approach allowed macroscopic and 
microscopic examination under changes of scale, shift and 
shape. The clustering and separation qualities are data 
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used, all 



dependent. Since no ship radar video data was 
observations were with respect to canonic classes. No 
strong groupings were detected. 



As alluded to earlier, the results for the exponential 
algorithm were less than satisfactory. Test shapes were 
carefully designed to minimize sampling effects, and to 
ensure low side lobes in the frequency domain. Many tests 
were conducted using each of the zero point correction 
methods with varied shapes, sample rates, and spectral 
resolutions. Classification on the basis of signal shape 
was very poor. The strongest correlation was between shapes 
of common duration. Table 1 shows a typical result of 
comparing a rectangular shape and a raised ramp. The 
scaling in each case was by 2 (60 degrees). Equation (36), 
Z=h(0)/ui was used, but the others offerred little 
improvement. Consistently, the strongest similarity was 
shown between shapes that had tbe same sample length, vice 
shape . 
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TABLE 1 



Canonic Shape Fourier - Mellin 
Feature Comparisons 



a. Peak Correlation Values 
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III. DIRECT MELLIN TRANSFORMS 



The last chapter developed a method of obtaining the 
Mellin transform by exponentially warping the signal prior 
to using an FFT block. This technique is referred to as the 
fast Mellin transform (FMT). Although the promise of scale 
invariant features is attractive, some of the problems 
encountered that make the FMT unattractive are reviewed 
here. The required sample rate varies with respect to the 
data, making general applications difficult. The tendency 
is to use more samples than required, which quickly becomes 
costly in an exponential sampling scheme. The need to 
exponentially warp (to interpolate) a set of new samples, 
is expensive in time required. For true scale invariance, a 
correction factor is required but because of the integral's 
unbounded nature at zero this correction factor is 
indeterminant. Several correction methods have been 
employed, but they depend on unspecified data 
characteristics. These effects combine to make the actual 
performance of the algorithm poor. Although scaling effects 
are mitigated, they remain an artifact, which is disturbing 
to classification attempts. This chapter outlines the 
effort to remove these limitations. Some useful Mellin 
properties are developed, and then applied to establish 
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several Direct Mellin Transforms (EMTs) which were built, 
and their performance compared. 

A. SONE USEFUL PROPERTIES 

Some general observations are made here about the Mellln 
transforms, and are followed by some specific relationships 
which were derived and applied. A property of the Mellin 
s-domain is that it is unaffected by scaling changes in the 
original x-domain . Figure 11a is a test shape in the 
x-domain. Two features are identified according to their 
amplitudes A and B, at x=a and x=b respectively. The ratio 
of a/b equals c. To be simply scaled by k, h(x/k) must 
remain the same in all aspects except that the distance 
between features has been changed according to the scale 
factor k. Figure lib shows the scaled domain. Features A 
and B are again identified at their scaled positions ka and 
kb. The signal property maintained by scaling is relative 
positional integrity, that is, ka/kb equals c as before. 
The positional integrity of the features, the ratio of 
their distance from the origin, to that of another's is 
unchanged. Restated, an operation 0[h(x)] in the x-domain, 
will leave the s-domain modulus invariant to simple scaling 
by k if 

0[A 0/A)J= (44) 
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Pure Scaling 
Figure 11 



Note that the entire domain is scaled, so no unsealed shift 
in the domain can he permitted as already discussed. In the 
liellin s-domain, any simple scaling in the x-domain results 
in a phase distortion (the modulus is invariant to k). 
Manipulations in the s-domain will leave that domain 
invariant to k, as long as the modulus is modified hy a 
multiplicative factor of constant phase. That is, if 'C(s)' 
is an arbitrary function of s (except that it does not 
depend on k) and i M [h ( x/k)] J= J H( s) j is also invariant to k, 
then their product is invariant to k and the x-domain 
remains simply scaled. For instance, in the x-domain the 
operator 0[h(x)3=x(h(x) ) , does not meet condition (44). 

7* -f (45) 

So the Mellin transform's modulus of xh(x/k) cannot be 
invariant to k. 

In Chapter II, a error was apparent because h(x) was not 
equal to zero for x=0. If h(x) could be modified with an 
operation that met the condition of equation (44) and 
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always produced a series that was zero at x=0 a general 
approach can he developed. Consider two operators, 



c> [!(«/*.)] = *r ( «■*<■*/*»* -c<.*/k) 

(46) 



cl 



o l Ac* /A)] * * 4^ (A C'X/A)) - -Pc 'x/A ) 



(47) 



Equation (46) will produce an acceptable f(i) as df/dx=0 at 
x=0. Equation (47) will always produce the required 
condition. To see the frequency domain equivalents, we must 
assume that the Mellin integral exists. Integrating by 
parts , 

I | * | Hi*) I 




jicx/Jk) 

cA 



4 sx 



j[ x s ~A.('*/a)J -s J A («/* ) * 



~l~ s j Ac*/a) I s (48) 

The result in the frequency domain is consistant with the 
conditions stated above. The limit of h(x) as x goes to 
zero or to infinity must be zero if equation (48) is to be 
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true. Similarly, under the same conditions, it can be shown 
that 



Equations (48) and (49) clearly do not apply where h(0) is 
not equal to zero. However, by using the x(d( h(x) )/dx) 
irodifying operation of (47), a series can always be zero at 
i=0. The function h(x) is further constrained by the fact 
that it must fall off faster than 1/x. This assumption must 
be valid and the modifying operator must be applied in the 
x-domain for the Mellin integral to exist in general. A 
Mellin transform of a function, after having a modifier 
applied, will be called a modified Mellin transform. 



The integral is close to a form which is realizable, except 
for the upper limit. For a finite sampled series, h(n) will 
be assumed zero outside of the interval 0^n<N. This 
truncation effects the transform of an otherwise infinite 
series. In this application the Mellin is applied to an FIT 
frequency spectrum. The truncation in the frequency domain 
is due to band limiting the signal prior to sampling. 
Scaling in the time domain will not result in simple 
scaling, but will add a dependence on the scale factor k 
that can not be removed by the transform. An approach by 




(49) 




(50) 
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Prost and Goutte is used to predict the size of the error 
[24,25]. First a suitable function will be selected and a 
relative error of truncation (RET) determined and applied 
to two scalings. The relative difference of the feature 
space is found and identified as the error. Remember ing 
that this is applied to a frequency spectrum, dh(f)/df is 
approximated by a function of the form. 



The modified Mellin transform of (51) over a finite 
frequency range would be approximately. 



The lower limit has been set in a manner to be consistant 
with Plancherel's theorem. A convenient worse case 
assumption is that the lower limit is essentially zero, but 
this depends, in general, on F and the data itself. 
Ha*(s,f) converges toward Ha(s), equation (50), in the mean 
square as F tends to infinity. The mean square error will 
not provide an expression for the error that can be used to 
correct for it but is useful to measure the effect of the 
truncation in more general terms. A relative error of 
truncation (RET) can be computed if the error is assumed to 
be distributed evenly over the range from which Ha*(s,f) is 




(51) 




(52) 
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computed. By employing Plancherel 's theorem, the mean 
squared values are, 




(53) 




(54) 



The RET is defined and solved as. 




(zf\F * t) * 



(55) 



But the limit F depends not only on the pass band F, but on 
the scaling in the f-domain. A relative error (RE) between 
a truncated spectrum and a scaled and truncated signal 
would be more complex, but worth the effort. 



where 0<k<l. Two observations should be made here. First, 
the relative error of truncation (55) and the relative 
error or difference between two truncated spectrums scaled 
differently (56) both depend on F. F depends on the cut off 
frequency of the low pass filter prior to any sampling. It 
is often chosen to minimize aliasing depending upon the 
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limitations of the sampling circuit. Second, RE depends on 
k as well. The importance of scaling differences to this 
data type can he readily seen. If equations (55) and (56) 
are valid, and if the range of k can be bounded (a design 
specification) F can be chosen to realize a stated RE. Or 
if F is fixed and k bounded, the RI may be determined for 
evaluation. If the data type is not appropriate, a more 
representative function may be determined and used in place 
of equation (51) to attain better expressions for RET and 
RE. 



3. ALGORITHM DEVELOPMENT 
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1. First Difference Approximations 



Although developed through a different rationale, 
this first algorithm was developed by Zwicke and Kiss [11] . 
Starting with a sampled series hi, the series is operated 
on by the modifier defined in equation (51), using the 
first backward difference to approximate the derivative 
with respect to x. Unit step size is assumed. 



Taking the trapezoidal rule to evaluate the modified Mellin 
integral (52) while recalling that h(0)=0, and h(N) is 
assumed zero, 



They can be calculated off line and stored to produce just 
the desired characteristics. The factor (2^m/M) could be 
any number that produces an interesting feature. Undesired 
features need not be computed (what in general was a N by M 
matrix, where M is the number of Mellin transform 
coefficients and N the number of x sample points). If a 
relatively small number of features is required, perhaps 




(57) 




(58) 



where s=-j2*Mn/M. The complex coefficients are 



CoS 2 ft' /m ) ~ £ J7 *\ /Art ) Jh*. C/*) 
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the processing will be manageable. Notice that no zero 
point correction is required. Cnly data changes contribute 
to the transform. These observations are valid for any of 
the modified Direct Mellin Transforms developed in this 
section . 

By using a central difference instead of the 
backward difference, a similiar result is obtained. 



Other numerical integrations may be used with improved 
results, and other methods can be used to increase the 
order of the approximation. 

To test the algorithms, a ramp and inverse ramp were 
used in a scaled and unsealed mode. The ramps and their 
scalings are shown in figure 13. Figure 14 is the analytic 
results of both waveforms plotted with the transform found 
using equation (60). This is a dramatic improvement over 
anything used with the methods discussed in the previous 
chapter. The noise of the signal appears to be diverging as 
frequency grows, but over the range plotted, the appearance 
is that of an algorithm trying to do well. 
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( 60 ) 
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2 . Second Difference Approximations 



A second difference algorithm can he achieved hy 
proceeding as before. The pure scaling operator used to 
prepare the signal is, 



Other second order operators have been used, and 
partitioned into forward and backward difference variations 
to (61), but this appears to be a basic and useful form. 
The color l/(s+l) is present to approximate the modified 
Mellin of equation (60). The term (s+1) is valid assuming 
that dh(x)/dx is exclusively upper bounded by ln(x)/x as it 
approaches zero or infinity. For comparison, another second 
difference algorithm was developed based on the modifier 
(x (d/dx ) (x (d/dx )h ( x) . 



This is roughly the sum of the methods defined in equations 
(60) and (62) above. The assumption for deriving the color 
1/s is even less restrictive than before, but the algorithm 
performs poorly compared to (62). These transforms depend 




And the new algorithm is. 




( 62 ) 







(63) 
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on second difference characteristics, but are modified by a 
1/s term which has a stabilizing effect. This is equivalent 
to a division by x and integration in the x-domain . The 
order of the approximation has been increased by the 
modifier. Results using equation (60) and (62) should be 
alike. Figure 15 is the results of using (62) compared to 
the closed form solution to the figure 13 test shape 
transforms. An improved performance over (60) is seen over 
some of the range, but a drop as frequency increases 
degrades the accuracy in figure 15b. 

3 . Higher Ilfference Approximations 

Higher difference approximations can be developed. 
For instance, one algorithm depending upon the third 
difference is. 



The performance of higher order algorithms becomes 
increasingly suspect because of the extreme weighting they 
apply to different parts of the series. Different 
algorithms exist, but this weighting is always a factor. A 
smoother transform is achieved, but a large error is likely 
to develop due to the algorithm's dependence on higher 
order derivatives and the nature of sampled data. However, 
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these comments are speculative since they were not 
confirmed experimentally. 



4 . Higher Order Integrations 

Higher order integration rules should he able to he 
used with a corresponding improvement in performance. One 
using the first difference with Simpson's rule was 
implemented with dissappointing results. Figure 16 is the 
result of such an implementation. The droop for the ramp 
input is apparent even though not present in the 
trapezoidal rule used in subsection 1 above. The higher 
frequency error is also more prevalent than before. A 
program error is of course suspected, but was never found. 
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IV. CLASSIFICATION PREPROCESS INC 



The problem of determining important signal 
characterises can be approached in at least two different 
ways. First, by trying to learn what is important in human 
recognition, and then trying to adapt a machine to emulate 
that behavior. Or second, by using successive transforms to 



remove information 


known 


to be 


superfluous 


to 


classification, while 


keeping 


enough 


informa tion 


to 



reliably assign an object to a class. Addressing the 
former, even though it is difficult to determine specific 
details, some key aspects of human visual recognition are 
discernable. Chief among these is that the intensity level 
of a scene, or object, does not appear to be as important 
as the relative position of the edges, i.e. the shape 
separating different intensities and frequencies [25-28]. 
Examples in scene analysis show clearly that the edges or 
shapes are far more critical to human recognition than the 
relative power differences themselves. The invariant shapes 
or angles in scene analysis find their analog in ratios 
between similiar points in differently scaled time series. 
Examination of many range only radar video ship signatures 
has provided a basis for noting that the relative position 
of time domain features remain constant over changes in 
aspect angle, while the relative intensity of the features 
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vary greatly as shown in figure 17. As the aspect changes, 
if this set of ratios remains constant within a ship class, 
then the set is identified as information worth preserving 
for the classifier. Conversely, since relative intensitiy 
is not a stable measure over aspect angle, or from among 
different ships of one class, that information should be 
intentionally removed to provide tighter natural 
clustering, with the minimum number cf features. 

A. INFORMATION REQUIRE! TO CLASSIFY 

There are two preconditions that must both exist for a 
set of possible input signatures to separate into distinct 
classes. First, the features of a particular class must 
have some common characteristic about them, and second, 
this characteristic must in some way be unique with respect 
to other classes. The assumption in existing radar 
signature classification projects is that there is enough 
Information in the signatures to permit this 
classification. Short of actually trying to classify with a 
set of realistic signatures, the analytical determination 
that sufficient information is present in a set of all 
possible ship signatures is difficult to approach. 

To establish how well the autonomous classifier is 
performing some measure of the classi f iabili ty of the set 



73 




A T I M E = 1 SEC 
ft B S < M A X > * 3 2 4 2 3 




A T I M E = 1 SEC 
RBS<f1flX>= 3 2766 

Ship Signatures 10 Degrees Apart 
Figure 17 



74 



received signals is desired. Tailing this, some discussion 
of the information capacity of the preprocessor should at 
least he considered. The preprocessor produces the TFT 
magnitude of a sampled signal as the output of the first 
stage. Most of the unique positional relationships of the 
signal upon which human recognition apparently depend has 
been destroyed. Next, the Mellin transfrom stage 
effectively distorts the signal and uses the magnitude of a 
second Fourier transform as the output features. Signals 
reconstructed on the basis of FFT phase information alone 
usually provide sufficient similarity to be associated with 
the original signal, whereas reconstruction on the basis of 
magnitude does not retain any significant detail except 
when the signal is symmetric [29] . The data is also masked. 
For most applications, the data is frequency shifted, 
filtered and sampled as a baseband signal. This windowing 
masks the magnitude characteristic and is specifically 
designed into the processing. The resulting FM features are 
insensitive to positional and scaling relationships that 
are necessary in human visual recognition. 



Analytic support is also available to quantify the 
importance of relative position of events [29]. By 

considering rms error (due to spectral phase and amplitude 
quantization for random signals) it has been concluded that 
approximately two more bits are required for the 



75 



quantization of phase information than amplitude for the 
same rms error. A separate analysis applied distortion rate 
theory to real-part, imaginary- part, and magnitude-phase 
encoding of the DFT of random sequences. The result was 
that phase required 1.4 hits more storage than magnitude 
for a similar error [30]. A third approach concluded that 
the Fourier phase includes 1.8 bits more information than 
the magnitude [31] . This was based on analysis of image 
reconstruction from kinoforms (phase-only holograms). The 
fact that phase-only reconstruction preserves much of the 
correlation between signals would suggest that the location 
of events tends to be preserved. Further, it seems that 
this information is lost by taking only the magnitude of 
the Fourier transform. Another interesting, albeit 
informal, view is apparent as one considers the phase-only 
signal as a spectral whitening process. 

F(f) = I F<f >!*->*' (65 

For reconstruction by phase alone, where the magnitude is 
set to one; 

( 66 ) 

Since the received radar signature will have an abundance 
of low frequency spectral lines and smaller high frequency 
components, the low frequency information is not weighted 
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as heavily as the higher frequency information in the phase 
reconstructed signal. This seems like it would accentuate 
sharp changes in the reconstructed signature without 
removing relative location information. The result is the 
summation of the different frequency components, all with 
zero phase (i.e., no positional or amplitude distribution 
information can remain). 



Considering the Kellin transform next, in a continuous 
case, the exponential warp does not lose any information. 
To zero the first data sample as required by the Chapter 
III modifications, will surely destroy information, but 
this may be confined to the EC term alone. The information 
lost during the final transform and magnitude is difficult 
to assess. Increased masking occurs due to the spectral 
truncation, so actual information loss may not be as great, 
but masking distortion may be greater than before. So 
approximately two bits of information are lost. Only a 
quarter of what was, remains. Information is also lost when 
the transforms are normalized in the processing so that any 
power calculation is also meaningless. Some interesting 
questions arise. After removing positional relationships, 
scaling, and power, what signal qualities remain and are 
they useful in classification? Although it is true that 
this insensitivity may add a certain robustness to the 
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system, the arbitrary loss of valid class if icat ion 
information should be minimized. 

By examining the quality that must be ignored, and by 
comparing its removal to what is actually removed by the 
processing, an interesting result will develop. The effect 
of a time domain shift on the frequency domain is an 
additive phase term, linearly related to the frequncy of 
the coefficient, as in equation (2). Most of the structure 
of the signal is held in the phase relationships with 
respect to the fundamental and higher frequency terms. So 
shift can be defined as the phase of the fundamental 
complex coefficient. 3y setting the phase to zero, and 
adjusting the other coefficients according to their 
component frequency, the structure of the signal is not 
lost but reconstructed about the fundamental as before. If 
there are N/2 spectral phase angles, only the fundamental 
needs to be zeroed to remove the shift. If the information 
is contained uniformly in the spectral phase, then only 2/N 
of this structural information needs to be removed. When 
the magnitude is taken to produce shift invarient features, 
all of the phase relationships are destroyed. (N-2)/N of 
the information once held in the phase was removed 
needlessly. The amount of information lost removing the 
shift can be made arbitrarily small. As N grows unbounded, 
the amount of information that needs to be removed tends to 
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zero. This result seems to be supported by experience. With 
enough samples of a shape, its position vith respect to the 
observation field is immaterial. In theory the principle is 
sound, but some practical limitations may degrade predicted 
performance. Recalling that an exponential warp translates 
scaling to shifting generalizes the result a bit further. 

4(t/A) = JL ( 67 ) 

The same priciple that permits simple shift removal is also 
valid for the removal of scaling dependence as well. Using 
the Beilin transform, scaling dependence may be removed by 
zeroing the fundamental and adjusting all the other 
coefficients as described above. For the FM preprocessor, 
the information lost removing the scaling and shifting 
dependence may be made arbitrarily small by increasing the 
number of spectral samples used. Since the number of 
spectral samples can be increased by filling zeros onto the 
finite signal in the original domain, this process does not 
effect the data sample rate. This approach was not verified 
experimentally, but represents a potentially powerful tool 
to analyse and improve the extracted feature space. 
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E. RANGE CNIY RADAR 



This section addresses ship classification on the basis 
of information gathered from a range only radar video ship 
signature. An example of such a signature has already been 
considered as figure 16. Classification by range only radar 
signatures is subject to the same distortions discussed 
above. Typically, the radar return is detected and isolated 
in a range gate that is sampled and digitized. The 
rectangular sampling window can be considered the range 
gate itself. The range gate is designed to ensure that the 
included range is greater than the maximum ship length so 
that the time/range windowing has no effect on the 
frequency spectrum of the signature, other than increased 
spectral resolution. The placement of the ship signature in 
the window is not set, neither from encounter to encounter, 
nor from pulse to pulse (jitter). The total effect joins 
together to produce the framing distortions. Signature 
scaling results from viewing the ship from different aspect 
angles. The sampling rate must be done at more that twice 
the inverse of the resolution of the receiver. Quantization 
levels are chosen in a manner to reduce that predictable 
random noise to an acceptable level. The pulse to pulse 
jitter is an ever present characteristic of the radar 
problem, but at a normal resolution (greater than 25 feet) 
integrating the return in the time domain removes the 
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jitter effect, reduces scintillation, and improves the 
resolution of the signature. Although predetection 
(coherent) integration is more efficient, post detection 
(noncoherent) integration is more common because of the 
convenience of not having to preserve the radar frequency 
(RF) phase. For post detection integration of n pulses, the 
signal to noise ratio would be something less than n times 
the signal to noise ratio for one pulse [32]. More 
important to the recognition problem itself, for a stable 
system by the law of large numbers [33], fluctuation of the 
average value of the return will be overcome. That is, with 
the integration of n pulses the resolution (R) will become 
finer as 



with a probability of 1-E , where S squared is the variance 
of the signal from pulse to pulse. For very high 
resolution, the cost of making the signature stable with 
respect to the integrator becomes prohibitive. It has 
become convenient to integrate the spectrum of the 
signature because the jitter effects can be completely 
removed. Another limiting factor in integrating over a 
period of time is that the position and aspect cf the 
target are dynamic. They change with time. A conceptually 
attractive solution is a recursive filter which weights the 




( 68 ) 
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Integrated pulses such that the older they becoire, the less 
weight is accorded to them. 



If the course and speed of the target ship are known 
from measurement of the target track, it is possible to 
infer the aspect of the ship. The range profile can give 
some estimation of size. Unfortunately, the three 
dimensional change in aspect angle, commonly suffered by a 
ship presents more than just a video signature scaling 
change. The radar cross section of even individual 
structural components of the radar target changes with 
respect to the aspect angle. The composite effect is that 
ship signatures vary greatly with aspect angle. The radar 
is an electromagnetic sensor, reacting to energy reflected 
from the target. These reflections are a result of 
scatterers that are related in dimension to the wavelength 
of the illuminating energy. Because of the great difference 
in wavelength between light and microwaves, what can be 
"seen” by radar may be quite different than that seen by an 
eye. Also, when measuring size or any distances with a 
radar of high resolution (less than 50 feet), an error can 
be incurred since the extremities of the target are not 
always good scatterers. Echoes from the forward or stern 
portions of the target might be observed in the noise, 
especially for a relatively low power radar [32], After 
reviewing hundreds of signatures, it appears that major 
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features, such as overall ship length and dominant mast 
structure are frequently dlscernable, hut vary in relative 
amplitude. Resonance, shadowing (one reflector hiding 
another), multipath returns, the mapping of three 
dimensional aspect changes onto a one dimensional time 
series, and the amplitude and phase of component returns 
summing constructively or destructively to cause a 
scintillation of the composite target. Some of the 
variations caused by these conditions can be lessened by 
integration but major effects remain causing the signature 
to vary in shape and content with the aspect angle. For 
this reason, a class feature volume cannot be reduced to a 
single point, but will remain a hypervolume in the feature 
space even in the ideal case. Any selection of features 
should try to minimize this volume. Features should be 
selected that are relatively insensitive to known 
superfluous effects. 

C. CLASS DISCRIMINATION 

In the last chapter, the major concern was removing two 
sources of variance with no classifying value. Algorithms 
were developed and canonic shapes generated to verify the 
algorithms and demonstrate the invariance to scaling. In 
the same manner, this chapter has reviewed information loss 
and process masking. The question of whether sufficient 
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information is preseat to classify was raised. To answer 
this question some simple classes of canonic figures are 
defined, and put through the entire FM preprocessor 
documented in Appendix E. The DMT algorithm found to offer 
the best scale factor rejection in Chapter III was used to 
generate the final features. The algorithm chosen is based 
on a second difference modification and is defined in 
equation (62), 



After canonic tests were made, preprocessor performance on 
several ship signatures was recorded. Although this was 
premature in the logical test sequence, the results are of 
some interest. 

1 . Test Shapes and Results 

Four test shapes were used. Figure 18 shows the test 
shapes. All were scaled and shifted originally to test for 
algorithm verification and demonstrate scale invariance. In 
this series of tests they were left fixed and used in 
different combinations to try to detect shape presence in 
the Ft* feature space. The object is to differentiate 
between different canonic classes. Figure 19 shows a 
comparison of the shapes, rectangle and triangle. A test 
combining the rectangle and triangle was the subject of 



W-/ 





(62) 
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figure 20. Finally in figure 21 a rectangle with two 
triangles is shown in the FM feature space. It is clear 
from the plots that rrost of the scale variance has been 
removed. Just as important, some quality does remain that 
differentiates between the canonic shapes. A square in 
general can be differentiated from a triangle. A ’’ship” 
with a single mass of superstructure can be separated from 
one with two such masses. 

Although it's clear from the plots that there is a 
unique quality left in the feature space to allow the time 
domain shapes to be classified, some quantified measure of 
system performance is required. The magnitude of the 
feeature vectors have all been normalized with respect to 
the first coefficient, so in that region little 
discrimination can be expected. For higher Beilin 
frequencies, noise dominates. A region of coefficients, 
11-100 was chosen to classify the shapes by correlation. 
The results are included as Table 2. The improvement in 
performance over that shown by Table 1 in Chapter II is 
dramatic. The methods in that earlier test resulted in the 
observed length of an object being the distinguishing 
criterion for classification. Using methods supported in 
Chapter III, variations due to scaling and shifting of the 
original domain have been removed. The features now reflect 
the shape of the object in the time domain. An unusual 
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effect is noted with respect to the difference squared 
analysis in Table 2h. Although TRI2 is more closely 
correlated to TRI1 than RECT2 the squared error shows just 
the reverse. The results are encouraging. The preprocessor 
has greatly simplified the classification problem for the 
canonic shapes above. 

2. Ship Signatures and Results 

A single ship was used to make these preliminary 
tests for this thesis. Signatures were taken every ten 
degrees around a ship from zero to fifty degrees. The 
results are plotted and compared over twenty degree aspects 
in figures 22-23. The signatures are the result of very 
high resolution radar signature data that has been degraded 
and smoothed to a lower resolution with essentially no 
noise present. Recalling that the purpose of the 
preprocessor was to remove variance due to pure shifting 
and pure scaling, leaving enough information for 
classification, to the eye there seems to be little 
encouragement from these results. It is recalled that a 
goal of this preprocessor is to make the classifiers job 
easier by removing dependence on shifting and scaling of 
the original data. The information that remains depends on 
unspecified signal characteristics that here appear to 
useful in discriminating shape classes and possibly ship 
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classes. However no real conclusion can be drawn at this 



point because of the snail data base and the absence of an 
automatic classifier to generate an optimal feature space. 
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V. CONCLUS IONS 



The preprocessor design began by considering a generic 
classification system. A distinction was drawn between the 
classifier and the preprocessor. The preprocessor is 
problem specific. It assists in the classification by 
extracting a set of features for the classifier. The 
extracted feature space has enhanced natural clustering on 
the basis of shape by removing information that was 
extraneous for classification. Two useless characteristics 
were identified as shifting and scaling. The preprocessor 
was designed to remove dependence on these two 
characteristics by using the invariant properties of a 
Fourier transform followed in series by a Mellin transform. 
The resulting set of coefficients are Fourier-Mellin (FM) 
f ea tures . 



A. REVIEW 

In Chapter II a Mellin transform was developed using the 
conventional digital processing approach which 
exponentially warps the domain and then transforms the 
spectrum by an FET . This method is sometimes identified as 
a fast Mellin transform (FMT). The unbounded behavior of 
the exponentially warped frequency spectrum was shown to 
result in a function that could not be transformed. That 
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is, the warped function has no transform because the Mellin 
integral is indeterminate at the lower limit. Error 
correcting techniques cannot compensate for the effect 
because the error itself cannot be computed in general. The 
bound for the error was found and seen to be the envelope 
for existing error correction functions. 

Chapter III developed some useful properties of the 
Mellin transform that were used to modify the signal so 
that the pitfalls isolated in Chapter II were avoided. This 
was done by modifying the input to the Mellin transform to 
always be transformable. To simplify the implementation and 
to control the effects of sampling, a direct Mellin 
transform was used for the development of the modifiers. 
Several suitable modifiers were determined and tested with 
differently scaled inputs for which the closed form 
solution was known. The direct Mellin algorithm that 
produced features closest to the closed form solution was 
chosen for use in the preprocessor. With the modifications 
in place, the preprocessor was tested and shown to produce 
features that were invariant to shifting and scaling. It 
was also shown that the features retained enough 
information to classify canonic shapes. 

Chapter IV discussed what type of information is 
required for classification. Signal structure or shape was 
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identified as key inf ormat ion . A discussion of wnat 
information is required for classification and a means of 
keeping the signal structure intact throughout the 
preprocessor was advanced, hut not empirically verified. 

B. FUTURE WORK 

The design FM preprocessor does produce a feature space 
with enhanced clustering, but problems remain to be 
resolved before the full potential of the system can be 
realized. There are three extant conditions that detract 
from the performance of the implemented preprocessor. 
First, the preprocessor is not computationally efficient. 
Second, most of the signal structure that should be vital 
to classification is obviously lost. And third, a complete 
verification of performance has not been conducted. The 
current preprocessor produces an enhanced feature base for 
a classifier, but attention to these main weak points will 
greatly improve the applied techniques. 

1 . Efficient Processing 

A direct Mellin transform using N spectral samples 
to transform into M Mellin spectral coefficients requires M 
by N multiplications. If only a few coefficients are to be 
used then the number of multiplications may be small. Using 
the FMT requires about N(20+ln(N)) arithmetic operations 
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for the interpolation and the FFT. If twenty-five or more 
Mellin coefficients are required, the FMT is faster. The 
modifiers used to prepare the spectrum for the direct 
Beilin transform will also work for as FMT, but this should 
be demonstrated experimentally. Because the FMT directly 
weights the lower numbered samples it may produce more 
accurate results as well. 

2 . Conservation of Information 

Chapter IV discussed the type of information 
required for visual pattern recognition in humans. The 
preservation of this information should be a specified 
design goal for the preprocessor. The preprocessors built 
for this thesis removed much of these vital signal 
properies. A means was introduced to limit or control the 
loss of structural detail by zeroing the fundamental phase 
and adjusting each of the remaining complex coefficients to 
reconstruct phase relationships. This may be done in the 
frequency domain as described, or by a similar operation in 
the time domain, shifting the centroid to zero. In either 
case, to transform more information about the signal will 
effectively increase the sensitivity of the features to 
characteristics in the time domain. This increase in 
sensitivity needs measured to confirm the approach. It is 
also possible that continued selective zeroing of 
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coefficients may offer improved, performance or robustness 
to the system as a whole. 

3 . Verification 

Although the improved performance was demonstrated 
with respect to ealier FM digital preprocessing, a direct 
improvement factor needs to be established. Time domain 
correlation should be used as a measure of original signal 
classif lability . This measure needs to be compared to the 
FP domain correlation recorded as Table 2 in Chapter IV. 
Next, the preprocessor or an improved version will have to 



be married 


to a class 


ifier and realistic 


data used 


to 


evaluate its 


effect on 


the classification 


system. 


The 


preprocessor 


built was 


design to be used 


on-line . 


An 



on-line classifier needs to be built as well. 

The systematic evaluation of ship profiles using FM 
features is still a requirement. For several ships, FM 
features must be singled out and plotted together as a 
function of aspect angle. The purposes are to establish a 
range of aspect over which classification may be possible, 
to evaluate changes of structural content as discussed in 
section two above, and to determine the beam signature as a 
classification "node”. Although these purposes assist in 
the system design, the third is more of an operational 
necessity. All ship signatures will degrade to the beam 
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aspect "node" so that this hypervolume in the feature space 
is occupied in common. Therefore the beam condition must be 
detected and withheld from ever entering the classifier. 
Nor should beam signatures be used for classifier training. 
The beam "node" needs to be determined separate from the 
classifier. Initially this information can be plotted to 
examine feature behavior and confirm the methodology, 
automated methods will quickly follow. These analysis 



funct ions 


may 


never reside in 


the 


clas 


sification system 


itself. 


but 


must be a part of 


the 


tool s 


used to develop a 



workable classification system. 
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APPENDIX A 



******* t***j*v*J*S:»>rJ** ****** ******** *****:*»*¥*<: **:»*»*** 

* THIS PROGRAM IS A DIGITAL IMPLEMENTATION OP * 

* A FAST FOURIER TBANSFOBM FOLLOWED BY A MELLIN * 

* TBANSFOBM. THE DIGITAL IMPLEMENTATION OF THE * 

* BELLI N IS AS AN EXPONENTIALLY SAMPLED SPECTRUM * 

* THEN HON THBOOGH AN FFT AND COBBECTED FOB THE * 

* EBBOB BY ANY ONE OF SEVERAL CORRECTION * 

* SUBROUTINES. COBCTX. * 

***<i*3(^***j|t:p*:(t* It******** ********************* ****** 



RECORDED PLOTS / PLOT NUMBER USING BECORD CALLS 
-TIME FUNCTION (CONT) / 1 
-SAMPLED TIME FUNCTION / 1 
-F FT (CONT) / 2 

-EXPONENTIALLY SAMPLED FFT / 2 
-UNIFORMLY SAMPLED FFT / 3 
-DISCRETE MELLIN FEATURES / 4 



RECOBDED PLOTS / PLOT NUMBER USING BECANG CALLS 
-TIME FUNCTION (CONT) / 1 
-FFT UNIFORMLY SAMPLED / 2 (MAG & PHASE) 
-MELLIN FEATURES / 3 (MAG & PHASE) 



150 

200 



DIMENSIGN XBEAL (200) , XIMAG (200) ,XINT (200) ,T (64) , 

1 ERT (200) 

ISCALE = 31 
MU=5 
M=2**MU 
N 0=7 
N=2**NU 

CALL SAMP (XREAL. XIMAG, N) 

CALL BECORD (XBEAL, XIMAG , N , ISCALE) 

COMPUTE THE ACTUAL SAMPLE POINTS AS PROVIDED BY THE 
SAMPLED VIDEO. 



N0=5 

N=2**NU 

CALL SAMP (XBEAL, XIMAG. N) 

CALL BECORD (XBEAL, XIMAG , N, ISCALE) 

ISCALE = 31 

LU=7 

L=2**LU 

EC 200 1=1, L 

IF (I.GT.N) GO TO 150 

XINT (I) =XBEAL (I) 

PBT7l)=XIMAjG(I) 

GO TO 200 
XINT (I) =0.0 
ERT (I)=0. 0 
CONTINUE 



101 



310 



300 

C 

C 

c 

c 

c 

c 

c 

c 

c 



400 

C 

c 

c 

c 

c 

c 

c 

c 

c 

399 



401 



415 

410 

C 

C 

C 

C 

C 

c 

c 

c 

ccc 



CALL FFT ( XINT ,PRT , L ,LU) 

CALL FFT(XREAL,XIMAG,N,NU) 

CALL RECORD (XREAL, XIMAG ,N *1 SCALE) 
CALL RECORD (XINT,PBT, L, ISCALE) 

ISC ALE = 31 
DO 300 1=1, N 
XBEAL (I)=SQRT 
XIMAG (lj=0. 

CONTINUE 



(XBEAL (I) **2+XIMAG <I)**2) 



CALCULATE THE NEW SAMPLE POINTS AND INTEBPOLATE 
TO FIND THE EXPONENTIALLY SAMPLED VALUES. 

CALL NUPTS (XINT,M,N) 

AfTEH NEW INTERVALS CALCULATED AND STOBED TEMPORARILY 
IN VECTOR XINT, STORE FOB LATER PRINTING IN ? VECTOR. 

DO 400 1=1,11 
PET (I) =0. 0 
1 (I) = XINT (I) 

CONTINUE 

WITH THE NEW TIMES IN VECTOR XINT, AND THE CURRENT 
SPECTRUM SAMPLES IN VECTCR XBEAL, COMPUTE THE NEW 
EXPONENTIAL SAMPLES AND ENTER THESE INTO VECTOR 
XINT BY USING THE CHOSEN INTERPOLATION METHOD. 
FINALLY, RECORD THE FIRST SAMPLE TO BE USED LATER 
TO CORRECT THE MELLIN TRANSFORM FOR LOW FREQUENCY 
LOST IN THIS EXPONENTIALLY SAMPLED TECHNIQUE. 

F0=XREAL (1) 

CONTINUE 

CALL INTP2 (XREAL, N, XINT ,M) 

WRITE (2,401)M 
FORMAT (14) 

CALL SMAX (XINT, PRT, M, SCALE) 

DO 410 1=1, M 
XREAL (I) = XINT (I) 

XINT (I) =SCALE * XINT (I) 

XIMAG (I)=0 

WRITE (2,41511 (I) , XINT (I) 

FORMAT (?10.5,7X, FI 0.5) 

CONTINUE 

SUBMIT THESE EXPONENTIALLY SAMPLED VALUES TO THE 
FINAL FFT BLOCK. 

CALL FFT (XBEAL, XIMAG, M, MU) 

APPLY THE CORRECTION TERM, FIND THE MAGNITUDE 
OF THE NOW SCALE AND TIME INVARIENT FEATURES. 

CALL CORCT1 (XREAL, XIMAG, M,F0) 

CALL CORCT2 (XREAL, XIMAG, M,F0j 
CALL RECORD (XREAL, XIMAG, M, ISCALE) 

STOP 

END 
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noon 



c 

c 

c 

c 

c 

c 



102 



101 

100 



103 



****************************************** 

* SUBROUTINE FFT, GIVEN THE COMPLEX * 

* SAMPLES (2**NU=N) , HILL RETURN THE * 

* N COEFFICIENTS USING A DFT. * 

****************************************** 



SUBROUTINE FFT (XHEAL, XIMAG, N , NU) 

DIMENSION XBEAL (N) , XIMAG (N) 

N2=N/2 

NU1=«U- 1 

K=0 



DO 100 L=1,NU 
DO 101 1=1, N2 
P=IBITR (K/2**NU1.NU) 
ARG=6.2&3185*P/FLOAT(N) 

C=C0S {ARG) 

S=SIN JARG) 

K1=K+1 
K IN 2=K 1 +N 2 

TREAL=XREAL-(K 1 N2) *C+XIMAG (K1N2) *S 



TIMAG=XIMAG jK1N2j *C -XBEAL (k 1 N2^ *S 
XBEAL (K 1N2) =XBEAL (K 1) -TREAL 
XIMAG (K 1N2) =XIMAG f K 1) -TIM AG 
XBEAL JR 1) =XBEAL (K T) +TREA1 
K 1) =3 



XIMAG 
K=K+ 1 
K=K+N2 



=XIMAG (K1) +TIMAG 



IF^K.LT.N) GO TO 102 

NU1=NU1-1 
N2=N2/2 
DC 103 K=1,N 
I=IBITR (K-1.NU) +1 
IF (I. LE . K) GC TO 103 
TREAL=XREAL (K) 
TIMAG=XIMAG JK) 

XBEAL (K) =XREAL (I 
XIMAG lK) =XIMAG (I 
XBEAL (I)=TREAL 
XIMAG [l\ =TXMAG 
CONTINUE 
RETURN 
END 



FUNCTION IEITR 
FUNCTION IBITR (J, NU) 
J1=J 
IEITR=0 
DO 200 1=1, NU 
J2=J1/2 

IEITR=IBITR*2+ (J1-2*J2) 
200 J 1= J2 

RETURN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 

c 

c 



*************************** ****************************** 

* INTP2 IS A SECOND OEDEfi INTERPOLATION BASED ON A 

* PCLINOMI AL EQUATING TO FIRST AND SECOND DERIVATIVES 

* APPROXIMATED BY CENTRAL DIFFERENCES. THE INPUT VECTOR 

* <X> CONTAINS UNIFORMLY SAMPLED DATA WITH UNITY 

* BETWEEN CONSECUTIVE SAMPLES. VECTOR <Y> IS INPUT 

* WITH THE NEW SAMPLE TIMES AND OUTPUT WITH 

* THE NEW SAMPLES. THERE ARE N UNIFORM SAMPLES IN <X> 

* AND M WARPED SAMPLES IN <Y>. 
********************************************************* 

SUBROUTINE^INTP2 (X A N, Y, Mj 



10 

25 



30 

60 



C 

C 

C 



C 

C 

C 



C 

c 

C 



DIMENSION X3 (3) ,X (N) , Y (Mj 



CHOSE THE X SAMPLE TIME CLOSES TO THE WARPED 
TIME HELD IN <Y> 

DO 60 IY=1,M 
DT = 100 
DO 10 I X= 1 . N 
IXTIME = IX - 1 
CTABS = ABS (DT) 

DIST = FLOAT (IXTIME) - Y (IY) 

DIST = ABS (DIST) 

IF (DTABS .LT. DIST) GO TO 25 
DT = Y (IY) - IXTIME 
CONTINUE 

IXTIME = IXTIME - 1 
J 1 = -1 
DO 30 J=1 ,3 
J 1 = J ♦ IXTIME - 1 
X3(J) = X (J1) 

CGNTINUE 
HUE = Y (IY) 

Y (IY) = FCN (X3 , TIME , IXT IME) 

CGNTINUE 
RETURN 
END . 

FCN IS THE INTERPOLATION RULE. 

FUNCTION FCN(X,T,I) 

DIMENSION X(3) 

TX = FLOAT (I) 

COMPUTE THE COEFFICIENTS. 

A = (X (3) - 2.* X (2) ♦ Xfin / 2. 

E = |x]3) - xnn / 2.0 - 2.* A * TX 
C = X(2) - A * IX** 2. - B * TX 

COMPUTE THE INTERPOLATED VALUE. 

FCN = A * (T**2) + B * T + C 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 

c 

c 



101 

100 

40 



30 

35 

60 



******* 9 ****************************************** 

* SUBROUTINE INTERP USES A LAGRANGE THIRD 

* ORDER METHOD OR A SECOND ORDER POLYNOMIAL 

* TC COMPUTE THE INPUT SAMELE WAVEFORM. 

* 

* THE SUBROUTINE INPUTS ARE: 

* X-SAMPLED WAVEFORM 

* N-THE NUMBER OF X SAMPLES 

* Y-INPUT SAMPLE TIMES 

* -OUTPUT ITERPOLATED SAMPLES 

* M-THE NUMBER OF I SAMPLES 
************************************************** 



SUBROUTINE INTERP (X , N. Y , M) 

DIMENSION X4(4) ,X(N) ,Y (M) ,SHFTX (156) 

FIRST CHOOSE XTIME NEAREST EACH YTIME, 

THEN COMPUTE THE INTERPOLATED VALUE AT THAT PT. 



DO 100 1=1, N 
15=1+5 

IF (I.GE. (N-5) ) 
SHFTX (15) =X (I) 
GO TO 100 
15=1- (N-5) 
SHFTX (15) =X (I) 
CONTINUE 
DC 40 1=1, M 
CONTINUE 
DO 60 11=1, H 
YP= 100 

DC 10 IX= 1 ,.N 
IXT=IX- 1 
YPABS=AES 



GO TO 101 



DIST=ABS (IX.I-Y (IY) ) 

IF (YPABS.LT .DIoT) GO TO 25 

YP= Y (IY) -IXT 

IX4=IX+5 

CONTINUE 

IX4=IX4-2 

DO 30 J=1 ,4 

X4 (J) =SHFTX (IX4+J) 

CONTINUE 

IM1=IX4+1 

10=1X4+2 

IP1=IX4+3 

IP2=IX4+4 

Y (IY)=YLAGR(X4, YP) 

CONTINUE 

RETURN 

END 
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FUNCTION YLIN MAKES A LINEAR INTERPOLATION 



FUNCTION YLIN(X4,YP) 
DIMENSION £4(4) 

IF (YP . LT. 0 . ) GO TO 50 
YLIN=X4 (2) + YP* (X4 (3) -X4 (2) ) 
HETORN 

YLIN=X4 (2) -YP* (X4 ( 1) -X4 (2) ) 

RETURN 

END 



FUNCTION YLAGR COMPUTES THE LAGRANGE MULTI-PLIERS 
AND MAKES THE INTERPOLATION FOR A CHOSEN OFFSET 
FROM THE CENTRAL SAMPLE X4(2) 



C 

C 

C 

C 

C 

C 

C 

C 

C 



100 



FUNCTION YLAGR (X4,YF) 

DIMENSION X4 (4) 

CM1=-YP*4YP-T) * (YP-2)/6.0 
C0= (YP**2-1)* (YP-2) /2.0 
CP1=- YP* (YP+1) * (YP-2) /2 . 0 
CP2=YP* (YP**2-1) /6 . 0 

YLAGR=CM1*X4(1) +C0*X4 (2) +CP1*X4 (3) +CP2*X4 (4) 

RETURN 

END 

** * **4*4* ******* *************************** ******* ****** 

* SUBROUTINE NUPTS CALCULATES THE M EXPONENTIAL SAMPLE 

* POINTS FROM THE N UNIFORM SAMPLES OF THE EXISTING 

* SPECTRUM IN PREPARATION FOR AH INTERPOLATION. THIS 

* EXPONENTIALLY SAMPLED SET OF POINTS ARE STORED IN 

* THE INPOT VECTOR X. 
******************************************************** 



SUBROUTINE NUPTS (X,M,N) 
DIMENSION £(156) 
UN=FLCAT(N)/2.0 + 1.0 
£M= FLOAT (M) -1.0 
DELZ=LOG (UN)./EM 
DC 100 1=1, M 
SI=FLOAT (I) -1.0 
VALUE = SI*DELZ 
X (IV=EXP (VALUE) 

CONTINUE 

RETURN 

END 
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C 

C 

C 

C 

C 

c 



0 

0 



******************************************************** 

* THIS IS A STUB PROVIDING A TEST SKIN RETORN 

* TC SIMULATE A SKIN RETURN. NORMALLY THIS SUBROUTINE 

* HILL ACTUALLY SAMPLE A SKIN RETURN. 

******************************************************** 

SUBROUTINE SAMP (XR, XI, N) 

DIMENSION XR(N) ,XI (N) 

DESIGNATE SCALE FACTOR AND UNSCALED TIME SHIFT. 

SCALE = 16. / 32. 

SHIFT = 0.0 / 32. 

SCALE = 1.0/SCALE 
T0= .5 + SHIFT 
N 1= N - 1 

EG 100 1=1, N 
XR (I) =0. 

XI(I)=0. 

CALCULATE THE TIME OF THE SAMPLE. 

SK= {FLOAT (I) -1.0) /FLOAT (N1) 

BUILD THE TEST SKIN RETURN. 

IN THIS (C 2-PER) CASE A DOUBLE PERIMID. 

TSCALE = (SK-T0) ♦SCALE 
IWFM = 1 



IF( 


IWFM .EQ. 1) 


GO 


TO 


10 




IF | 


IWFM .EG. 2) 


GO 


TO 


50 




W = 
IF 


8. / 32. 
(TSCALE .LT. 


-W) 


GO 


TC 


100 


IF 


TSCALE .GT. 


+ wj 


GO 


TO 


100 


IF 


(TSCALE .LT. 


TO) 


GO 


TC 


20 



XR(I) 
GO TO 



10 



i 



TO - TSCALE) * 10.0 



X R (I) = W * 10.0 - (TSCALE - TO) * 10.0 



GO TO 100 
CONTINUE 
THE FOLLOWING 



WAVE 



FORM IS A SQUARE WAVE 16 SAMPLES 
WIDE CENTERED AROUND THE SIXTEENTH SAMPLE. 

SCALING AND SHIFTING DONE ABOVE WILL EFFECT 
WAVEFORM ACCORDINGLY. 



THE 



TSCALE = (SK-TO) *SCALE 
EDGE = 1. / 4. 

IF (TSCALE .LT. -EDGE) GO TO 100 
IF (TSCALE .GT. +EDGE) GO TO 100 
XB (I) = 1.0 
100 CONTINUE 
RETURN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



[****************************************************** 



c 

c 



** 

* THIS SUBROUTINE RECORDS A SELECTED DATA SET, 

* ' - -- - - ~ 

* 

* 

* 

* 

* 



AND SCALES THE INDEX TO 32 SAMPLES (0-31). 
THE INPUTS ARE XREAL* XIMAG L AND N ]THE 
NUMBER OF SAMPLES) 

VARIAELE. 

IF ISCALE = 0 



'isCALE Is A CONTROL 



* 

* 

* 

* 

* 

♦ 

** 



50 

40 

60 



75 

100 



100 



. NO SCALING IS PERFORMED 
SCALE - 1, AXIS = 31 
SHIFT = 1 (STARTS AT T = 0 ) 

ISCALE < 0 SHIFTING AND SCALING OCCURS 

| ISCALE j = AXIS LEN. A SHIFT IS INCORP- 
ERATED SO THAT AXIS STARTS AT 1. 

ISCALE > 0 AXIS SCALING OCCORS BUT THERE IS 
NO SHIFT (IE., THE AXIS STARTS AT 0.) 

*********************************************** ******* 

SUBROUTINE RECORD (XR, XI, N, I SCALE) 

DIMENSION XH(N) , XI (N) ,REC (200) 

SCALE=1 .0 
AXIS = 31 
SHIFT = 1 
WRITE (2,50) N 
FORMAT (14) 

IF (ISCALE .EQ. 0) GO TO 40 
AXIS = FLOAT (ISCALE) 

AXIS = ABS (AXIS) 

IF ( ISCALE .LT. 0) SHIFT = 0 
CALL SMAXJXB, XI, N, SCALE) 

DO 100 1=1, N 

SI= (I-SHIFT) *AXIS/ (FLOAT (N) - 1.0) 

HEC jl) =SQE T(Xl(I) **2+XR (1)**2) 

CCNT INU E 

EEC (I) = SCALE*REC (I) 

WRITE (2,75) SI.REC (I) 

FORMAT (F10. 5, 7X, FI 0.5) 

CONTINUE 

RETURN 

END 



SUBROUTINE SMAX (XR , XI, N , SCALE) 
DIMENSION XB(N) , XI (N) ,T (200) 
XMAX= 1.0 
DO 100 I=T, N 

T (I).=SQBT 'XR(I) **2 + XI (I)**2) 

IF (T (I) .LT. XMAX) GO TO 100 



(T I) 

XMAX = T (I) 
SCALE=1 . O/XMAX 
CONTINUE 
RETURN 
END 
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c ********************************* ****** ****** ****** 



C * THIS CORRECTION SUBROUTINE USES ONE OF * 
C * TflE SIMPLER CORRECTIONS FOR THE MELLIN * 
C * TRANSFORM. THE CORRECTION IS A PURE IMAGINARY. * 
C * * 
C * CORRECTION = -FO/OMEGA * 
C * * 
C * THEN A MODIFICATION IS MADE TO THE ENTIRE * 
C * TRANSFORM BY A MULTIPLICATION BY OMEGA. * 



C *************************************************** 



SUBROUTINE CORCT1 (XR,XI,N,F0) 
DIMENSION ZB (N) ,XI (N) 

DO 100 1=1, N 

OMEGA = FLOAT (I) - 1. 

XR (I) = XR (I) * OMEGA 
IF (I .EQ. 1 ).GO TO 90 
XI (I) = CXI (I) - FQ/OHEGA) * OMEGA 
GC TO 10O 
90 XI (I) = FO 



100 CONTINUE 

RETURN 

END 



C 

C 

C 

C 

C 

C 

C 

C 



******************************************************* 

* THIS CORRECTION APPLIES THE MORE COMPLECATED 

* EXPRESSION, 

* CORRECTION = F0/2 + JCOT (FO/OMEGA) 

* AND THEN MODIFIES THE ENTIRE TRANSFORM BY 1/OMEGA. 

******************************************************* 



SUBROUTINE CORCT2 (XR, XI ,N ,F0) 

DIMENSION XR (N) , XI ( N) 

CO 100 1=1, N 

CMEGA = FLOAT (I) - 1. 

XR (I) = (XR- (I) + FO/2.) * OMEGA 

IF (I .EQ. 1 ) GO TO 90 

XI (I) = (XI (I) - (FO/2. ) ♦COTAN (OMEGA) ) IOMEGA 
GO TO 100 

90 XI(I) = FO/2. 

100 CONTINUE 

RETURN 

END 
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APPENDIX B 



************ ♦M***«^*W**#W#*»'(i***«****>li»a*?*?***4iW**» 

* THIS PROGRAM, FOURIER DIHECT MELLIN, TAKES AN INPUT 

* WAVEFORM FROM LOGICAL DEVICE 2, PERFORMS AN FFT 

* FOLLOWED BY A MDHT CUTPUTING THE FEATURES TO LOGICAL 

* DEVICE 3 FOR LATER PLOTTING BY MELPLT. THE MAJOR 

* DATA STRUCTURES AND SUBROUTINES ARE LISTED BELOW 

* IN ORDER OF THEIR APPEARANCE. THE SUBROUTINES ABE 

* DESCRIBED IN MORE DETAIL WHERE THEY ARE ACTUALLY 

* LISTED IN THE PROGRAM. 
********************************************************* 



MAJOB DATA STRUCTURES: 

< WFM> - THE INPUT WAVEFORM (REAL) 

<XFM> - THE MELLIN TRANSFORM (COMPLEX) 

<CPHI> - REAL MELLIN COEFFICIENTS 
<SPHI> - IMAGINARY MELLIN COEFFICIENTS 
<STAND> - AN ARRAY HELD FOR LATER COMPARISON 
OR OTHER USE. 

<PHT> - AN ARRAY NORMALLY USED TO HOLD REAL 
DATA TEMPORIL Y . A WORK SPACE. 

<IXT> - X AXIS TITLE FOR PLOTTING 
<IYT> - Y AXIS TITLE FOR PLOTTING 
<KEY> - NUMBER OF PLOT THIS GRAPH 



SUBROUTINES: 

WAVE - READS AN ARRAY FROM LOGICAL DEVICE 2, 
AND FILLS ZEROS TO MAKE A TOTAL OF 256 
SAMPLES. THE OUTPUT IS IN <SFM>. 

FFT - AN FFT BLOCK 

COEF - COMPUTES THE MELLIN TRANSFORM SAMPLE 
WEIGHTS. THESE ARE COMPLEX 
NUMBERS WHOSE REAL AND IMAGINARY 
PARTS ARE STORED IN <CPHI> AND 
<SPHI> RESPECTIVELY. 

DMTM - APPLIES A MODIFIED DIRECT MELLIN 
TRANSFORM TO AN INPUT WAVEFORM 
PUTTING THE OUTPUT IN <XFM>. 

THE ALGORITHM IS BASED ON A FIRST 
BACKWARD DIFFERENCE. 

SMT - APPLIES A MODIFIED MELLIN TRANSFORM 
BASED ON A SECOND DIFFERENCE. 

SMT2 - APPLIES A MODIFIED MELLIN TRANSFORM, 
DIFFERENT THAN SMT, BUT ALSO BASED 
ON THE SECOND DIFFERENCE. 

CDMT - APPLIES A MODIFIED MELLIN TRANSFORM 

JUST AS DMTM ABOVE, EXCEPT THAT THE 
CENTRAL DIFFERENCE IS USED. 
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SIMP - APPLIES A dO DIFIEO MELLIN TRANSFORN 
USING A BACKWARD DIFF ERZNCE AS IN 
DMTM. EXCEPT THAT THE INTEGRATION 
IS BY SIMPSON'S ROLE INSTEAD OF THE 
TRAPEZOIDAL ROLE. 

XA3 - TAKES THE MAGNITUDE OF THE COMPLEX 
TRANSFORM <XFM> AND POTS THE 
MANITODE IN A SPECIFIED VECTOR. 

STOW - NORMALIZES A VECTOR BY ITS 

MAGNITUDE AND WRITES IT TO 
LOGICAL DEVICE 3 WITH A TITLE 
FROM LOGICAL DEVICE «*. 

HOLD - LOADS ONE VECTOR INTO ANOTHER. 

ALTER - CHANGES <HFM> BY SCALE &/OR SHIFT 
AND OUTPUTS TO A SPECIFIED ARRAY 

INTP3 - A SECOND ORDER SPLINE INTERPOLATION. 

CFOHH - PROVIDES TWO CLOSED FORM SOLUTIONS 

FOR VERIFYING THE MELLIN ALGORITHMS. 

TITLE - ENTITLES THE PLOTS ON THE BASIS OF 
THE CALLING PROGRAM. 



Ill 



c 



c 

10 



c 

c 

c 

c 

c 

c 



50 



c 

c 



c 

c 

c 

100 



c 



c 

c 



200 



c 

c 

c 

c 



SO THE MAIN PROGRAM STARTS! 

DIMENSION PET <256) , WFM (256) , STAND (256) , IF (5) 
COMMON XFM (256.2) ,CPHI (256, 12b) ,SPHI (256, 126) , PI, 
1IXT(10)-,IIT(ief ,k£i 

DATA IF/* FR* , ' EQOE* , * NCI *,* *,* •/ 

PI = 3. 141592654 



HOB MANY WAVEFORMS ARE TO BE TRANSFORMED? 
BEAD (2, 10) NUMSFM 
FORMAT (14) 



STMS IS THE NUMBER OF TIME SAMPLES (INCLUDING 
ANY ZERO FILLING) . IT IS A POWER OF TWO 
FOR THE CONVENIENCE OF THE FFT. 

BPTS IS THE NUMBER OF SAMPLES INPUT TO THE 
MELLIN TRANSFORM BLOCK. THE COEFFICIENTS 
ARE COMPUTED NOW. 

NU = 8 

STMS = 2**NU 
MPTS = NTMS/2 
FORMAT (14) 

NCOEF = NTMS 

CALL COEF (NCOEF, MPTS) 

SET UP THE LOOP FOR THE NUMBER OF WAVEFORMS 
TO EE PROCESSED. 

DO 500 IWAVE=1,NUMWFM 

GET THE NEXT INPUT WAVEFORM. 

CALL WAVE (WFM, NTMS) 

ZERO THE <STAND> VECTOR TO BE USED AS THE 
IMAGINARY PART OF THE NTMS TIME SAMPLES. 

DO 100 1=1, NTMS 
STAND ay = Q.O 
CONTINUE 

CALL STOW (WFM, NTMS) 

TAKE THE FFT. 

CALL TITLE (IF) 

CALL FFT (WFM, STAND, NTMS, NC) 

TAKE THE MAGNITUDE AND PUT IT INTO THE COMMON 
WAVEFORM <WFM> . 

DO 200 1=1, NTMS 



JJVJ X— * f lUUJ 

WFM (I) = WFM (I) ** 2 + STAND(I)**2 
”~Mjl) = SQBT (WFM (I)) 

NTINUE 



WFM 
CON 

CALL STOW (WFM, NTMS) 

TAKE THE MEiLLIN TRANSFORM OF THIS SPECTRUM 
USING THE FIRST HALF OF THE FFT SAMPLES , 
OTHERWISE KNOWN AS MPTS=NTMS/2. THESE ARE THE 
ONLY UNIQUE VALUES. 

CALL DMTM (WFM, MPTS, NCOEF) 

CALL XAB (PRT, NCOEF) 

CALL STOW (PRT, NCOEF) 



112 



onnoo nnnnfi oo nnnnnnnrm 



NEXT CALL THE SECOND ORDER ROLE DMT 
SOBROOTINES. BOTH COMPOTE THE MELLIN 
OSING THE SECOND DIFFERENCE APPROXIMATION 
INSTEAD OF THE FIRST DIFFERENCE APPROXIMATION 
ABOVE. OTHERWISE THE APPROACH IS THE SAME. 



CALL SMT (WFM,MPTS, NCOEF) 

CALL XAE(PRT,NCOEF)_ 

CALL STCW (PRT, NCOEF) 

CALL SMT2 (WFM,MPTS, NCOEF) 

CALL XAB (PR-T, NCOEFf 
CALL STCW (PBT, NCOEF) 

CALL CDHT WHICH USES THE CENTRAL DIIFERENCE 
ROLE FOR APPROXIMATING THE TRANSFORM. 

CALL CDMT (WFM,MPTS, NCOEF) 

CALL XAB (PRT.NCOEFi 
CALL STCW (PRT, NCOEF) 

CALL SIMP WHICH COMPUTES THE MELLIN TRANSFORM 
OSING THE FIRST DIFFERENCE ALGORITHM AND 
SIMPSON'S ROLE TO COMPOTE THE MODIFIED 
MELLIN TRANSFORM. 



CALL SIMP (WFM,MPTS, NCOEF) 

CALL XAB (PRT, NCOEF) 

CALL STCW (PRT, NCOEF) 

END THE TRANSFORM LOOP. THE TRANSFORM HAS BEEN 
COTPOT TO LOGICAL DEVICE 3 AND PREPARED WITH 
TITLE INFORMATION PROVIDED BY LOGICAL DEVICE 2 
FOR PLOTTING WITH PROGRAM MELPLT FORTRAN. 

STAY IN THE LOOP IF MORE WAVEFORMS ARE AVAILABLE. 
500 CONTINUE 

CALL CFORM (PRT, NCOEF) 

CALL CFORM (ERT, NCOEF) 

STOP 

END 
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*****************************.?********************* 

* THE WAVE SUBROUTINE PRODUCES A WAVEFORM, IN 

* BEAD FROM THE LOGICAL UNIT 2, NORMALLY A SHIP 

* DATA FILE. THE N SAMPLES ARE READ AND THEN 

* ZEROS ABE STUFFED TC Fill THE 256 SAMPLES. 

* THE REQUESTED WAVEFORM IS OUTPUT IN THE 

* COMMON ARRAY WFM (256) . 



SUBROUTINE WAVE (WFM , NPT S) 

DIMENSION WFM (NETS) , ID (5[ 

COMMON XFM (25(3 , 2) ,CPHI (256,128) ,SPHI (256, 128) ,PI, 
1IXT (10) ,IYT (10) ,KEY 

DATA ID/' RAN ' , ' GE S' , ' AMPL • , * ES •/ 

READ (2, 10) KEY 



10 FCRMAT (14) 

20 FORMAT ?F1 0.5) 

IF ( N .EQ. NETS) RETURN 
N = N + 1 
DO 100 I=N, NPTS 
WFM (I) =0.0 
100 CONTINUE 
RETURN 
END 



30 

50 
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******************************:*»* ********* 

* SUBROUTINE EFT , GIVES THE COMPLEX * 

* SAMPLES ( 2**N0=N) , WILL BETUEN THE * 

* U COEFFICIENTS USING A EFT. * 

************* ****** *********************** 



02 



01 



00 



03 



SUBBCUTINE FFT (XB EAI, XIMAG, N r NO) 
DIMENSION XBEAL (N) , XIMAG (N) 

N2=N/2 
N 01 = NO- 1 
K=0 

DO ICO L=1,NU 
CO 101 1=1. N2 
P=IBITR JK/2**NU1 , NO) 
ABG=6.28318'5*P/FLOAT(N) 

C=COS (ABG) 

S=SIN <ABG) 

K1=K+1 
K 1N2=K1 +N2 

TREAL=XBEAL (K1N2) *C+XIMAG (K1N2) *S 
TIM AG=XIMAG (K1 N2) *C-XBE AL (l 



K >N2) *S 

XBEAL (K 1N2) =XREAL (KI)-TREAL 
K1N2) =XIMAG fKll-TIMAG 
K 1) =XBEAL (K1) ♦THE AL 
K 1) =XIM AG (K 1) +TIMAG 



GO TO 102 



,N 



XIMAG 
XBEAL 
XIMAG 
K=K+ 1 
K=K+N2 
IFjK.LT.N) 

N01=N01— 1 
N2=N2/2 
DO 103 K= 1 , 

I=IBITR (K-1.NO) +1 
IF/I.LE.K)GC TO 103 
TREAL=XBEAL(K) 
TIMAG=XIMAG iKl 
XBEAL (K)=XBEAL (_ 
XIMAG (K)=XIMAG (I 
XBEAL (I j =TBEAL 
XIMAG (I) =TIMAG 
CONTINUE 
BETUBN 
END 



,; i 



FONCTION IBITR 
FUNCTION IBITR (J, NU) 
J1=J 
IBITR=0 
DO 200 1=1, NO 
J2=J 1/2 

IEITR=IBITB*2+ (J1-2*J2) 
200 J1=J2 

BETORN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



************************************************** 

* THE COEF SUEROUTINE CONFUTES THE MELLIN 
COEFFICIENTS IN TO COMMON ARRAYS CPHI AND 
SPHI THAT REPRESENT THE REAL AND IMAGINARY 
FARTS RESPECTIVELY. THE TERMS ARE COMPUTED 
EY THE FORMULA: 



* 

* 

* 

* 

* 

**:M*** ****** ************************************* 



PHI (I/O) = J**S , WHERE S = A NORMALIZED 

DISCRETE RADIAN FREQUENCY 



200 

100 



SUBROUTINE COEF (NCOEF. NETS) 

COMMON IF M (256 f 2) ,CPHI<256, 128) ,SPHI (256, 128) ,PI, 
1IIT (10) ,IYl (10) ,KEY 
DO 100 I = 1, NCOEF 
HI = FLOAT (II 

CMSGA = 2.0 * PI * HI / 36. 

DO 20C J = 1, NPTS 
RJ = ELOAT (J) 

CPHI (I, J) * COS (OMEGA * ALOG (RJ) ) 

SPHI (I. J) = SIN (OMEGA * ALOG (RJ) ) 

CONTINUE 

CONTINUE 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



******************************************************* 

* THE DMTM SUBROUTINE PERFORMS A DESCRETE HELLIN 
TRANSFORM ON THE ARRAY BFM. THE FORMULA 
FOR ONE MELLIN FREQUENCY VALUE IS: 

XFM (I) =SUM (K=1 TO NETS) (WFM (1+ 1 ) -WFM (I) ) *K**S 



THE COMPLEX COMPONENTS CF K**S ARE COMPUTED PRIOR 
TC CALLING DMTM AND STORED IN THE COMMON ARRAYS 
CPHI AND SPHI AS REAL AND IMAGINARY PARTS 
RESPECTIVELY. THE ALGORITHM IS BASED ON THE 
FIRST DIFFERENCE. THE TRAPEZOIDAL RULE IS USED 
FOR THE INTEGRATION. THE COMPLEX OUTPUT FOR 
THE TRANSFORM IN IN THE COMMON ARRAY <XFM>. 



******************************************************* 



SUBROUTINE DMTM (SAMP, NPTS,NCOEF) 

DIMENSION SAMP (NPTS) ,ID(5) 

COMMON XFM (256 , 2) ,CFHI (256,128) ,SPHI (256, 12 8) , PI, 
1IXT (lOJ^IYIJIO) jKEY 



DATA 



)/• D 



LLIN' , ' FRE* 



, 'QUEN* , ' CY 



/ 



CALL TITLE (ID) 

DC 100 1= 1 , NCOEF 
XFM (I,1)=0.0 
XFM fl,2j=0.0 
N 1 = NPTS - 1 
CO 200 J= 1 , N 1 
J1 = J + 1 





XFM 


1,1 


= CPHI 


(I,J) 


* (SAMP(j; 


) - SAMP 


( J 1) ) 


♦ XFM (1,1) 


200 


XFM 

CONI 


i, 2 i 

?INUS 


= SPHI 


ii,j> 


* (samp(j! 


| - SAMP 


U i)) 


♦ XFM (1,2) 



100 



CONTINUE 

RETURN 

END 
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C 

C 

C 

C 

C 

C 

c 

c 

c 

c 



***3r#****3M‘#**#***>j'**************:*#’ft$**#******4*:p4>****¥* 

* SMT IS A SUBROUTINE THAT COMPUTES A NUMERICAL 

* APPROXIMATION TO THE MELLIN TRANSFORM AS DOES 

* EMTM ABOVE. SMT USES A TRAPAZOI DAL APPROXIMATION 

* TC COMPUTE THE TRANSFORM, BUT USES THE SAME 

* COEFFICIENT MATRIXES <CPHI> AND <S?HI> CONTAINED 

* IN COMMON. 

***'4^*****4:M********3M‘$**************sft*****9:l'********** 



SUBROUTINE SMT (SAMP ,NPTS.NCOEF) 

DIMENSION SAMP (NPTS) ,ID j$) 

COMMON XFM (256,2) ,CPHI (256, 128) , SPHI (256, 128) ,PI, 
1IXT(10l,IYT(10) . K&Y 

CATA ID/* S-ME* , f LLIN' , * FRE ' , » Q DEN* , ' CY '/ 

CALL TITLE (ID) 



INITIALIZE IHE INPUT ARRAY AND COMPUTE 
THE LOOP CONSTANTS. 

N 1 - NPTS - 1 
N2 = N 1 - 1 

SET UP THE TRANSFORM LOOP. THE OUTER LOOP 
SETS UP THE COEFFIECIENTS WHILE THE INNER 
LOOP COMPUTES THE SUM WHICH ARE THE 
COEFFICIENTS. 



00 



DO 200 J 
XFM (J,1 
XFM (J ,2 
DC 100 

10 = I 

11 = I + 

12 = I + 



= 1, NCOEF 
0 
0 

N2 



in; 



1 

2 



DELTA = SAMP (10) - 2 .* 
XFM (J,1) = XFM (J, 1) 
XFMjJ,2i = XFM (J, 2) 
CONTINUE 



SAMP (II) + SAMP (12) 
+ DELTA*CPHI (0,1) *1 
♦ DELTA*SPHI (0,I)*I 



XFM (J,1) = XFM (J, 1) 
XFM (J,2) = XFM (J, 2 ) 



/ (FLOAT (J) *PI/18.) 
/ (FLOAT (J) *PI/18.) 



XFM ( J , 1 ) = XFM (0,1) / SQRT(1+ (FLOAT (J) *PI/18.)**2) 
XFM (0.2) = XFM (0,2) / SQRT ( 1+ (FLOAT (0) *PI/1 8. ) **2) 
200 CONTINUE 
EETURN 
END 
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noon no 



C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



4 * 4 * 4 * 4 4 * *** 4 *** 4 * 44 * 4 * 4 ***** 4 * 4 * 4 ** 4 * 4 * 444 * 4 * 444 * 44444 * 

* SMT2 IS A SUBROUTINE THAT COMPOTES A NUMERICAL 

* APPROXIMATION TO THE HELIIN TRANSFORM AS DOSS 

* EMTM ABOVE. SMT OSES A TRAPAZOIDAL APPROXIMATION 

* 1C COMPUTE THE TRANSFORM, BUT USES THE SAME 

* COEFFICIENT MATRIXES <CPHI> AND <SPHI> CONTAINED 

* IN COMMON. 

* 44 4 * ***************************************** ********** 



SUBROUTINE SMT2 (SAMP, NPTS , NCCEF) 

DIMENSION SAMP (NPTS) -ID (5) 

COMMON XFM (256, 2) ,CPHI (256, 128) , SPHI (256, 128) , PI, 
1IXT (10) ,IYT (10) -KEY 

DATA lb/’ S2-M 1 , * ELL I* , * N FR 1 , 'EQUE* NCY •/ 

CALL TITLE (ID) 



100 

C 

200 



INITIALIZE THE INPOT ARRAY 
THE LCOP CONSTANTS. 

N 1 = NPTS - 1 
N 2 = N 1 - 1 



AND COMPOTE 



SET OP THE TRANSFORM LOCP. THE OUTER LOOP 
SETS OP THE COEFFIECIENTS WHILE THE INNER 
LCOP COMPOIES THE SOM WHICH ARE THE 
COEFFICIENTS. 

DO 200 J = 1 , NCOEF 
XFM (J , 1 ) = 0.0 
XFM (0,2) = Q.O 
DO 100 I = .1, N2 

10 * I 

11 = I ♦ 1 

12 = I + 2 

DELTA = I* (SAMP (10) - 2.* SAMP (II) * 

DELTA = DELTA ♦ (SAMP (12) - SAMP (10) 

XFM (J , 1 ) = XFM (J, 1 ) * DELTA*CPHI (J 
XFMjj,2i = XFM (J, 2 ) ♦ DELTA*SPHI (J 



CONTlfiU] 

REMOVE THE COLOR 



•8 



SAMP (12) ) 



XFM(J,1) * XFH(J,1) / SQRT (1+ (FLOAT (J) * PI/ 1 8 . ) **2) 
XFMJJ.2J = XFM(J,2) / SCRT ( 1 + (FLOAT (J) *PI/18. ) **2) 
CONTINUE 



RETURN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



************************** * 4 ** ************************* 

* THE CCHT SUBROUTINE FERFORMS A DESCRETE MELLIN 
TRANSFORM ON THE COMMON ARRAY »FM. THE FORMULA 
FOR ONE MELLIN FREQUENCY VALUE IS: 

XFM (I) =SUM (K= 1 TO NPTS) (WFM (I ♦ 1 ) - «FM (I- 1 ) ) *K**S 



THE COMPLEX COMPONENTS OF K**S ARE COMPUTED PRIOR 
TO CALLING CDMT AND STORED IN THE COMMON ARRAYS 
CPHI AND SPHI AS REAL AND IMAGINARY PARTS 
RESPECTIVELY. THE CENTRAL DIFFERENCE IS USED. 



******************************************************* 



SUBROUTINE CDMT (H, N ETS , NCOEF) 

DIMENSION H (NPTS) , ID (5) 

COMMON XFM (256,2) ,CPHI (256, 128) ,SPHI (256, 128) ,PI, 
1IXT (10),IYT (10) -KEY 

DATA ID/* C-ME * , ' LLI N ' , * FRE * , ' QU EN‘ , • CY •/ 

CALL TITLE (ID) 



C 

c 



c 



100 

c 



200 



SET UP DERIVATIVE OF INPUT VECTOR <H(I)> 
IDENTIFIED HERE AS <G>. 

NG = NPTS - 1 

DO 200 J = 1 , NCOEF 
XFM (J,1) = 0. 0 
XFM(J,2) = 0.0 

OMEGA = FLOAT (J) * PI / 18. 

EC 100 I = 1 , NG 
IH1 = I 
IP1 =1+2 

DH = (H (IP1) - H (IM 1) ) / 2. 

COMPUTE THE J-TH COEFFICIENTS BY THE SOM. 

XFM (J , 1 ) = XFM (J, 1 ) + DH * CPHI (J,I) 

XFM (J,2) = XFM (J, 2) + DH * SPHI (J, I) 

CONTINUE 



COLOR = ‘ 
XFM (J,1) 
XFM (J , 2 ) 

CONTINUE 



BETURN 

END 



XFM (J, 1) * COLOR 
XFM (J, 2) * COLOR 



120 



# * * ****** 



c 

c 

c 

c 

c 

c 

c 



***************************************** ******** ****** 

* THIS SUBROUTINE USES SIMPSON'S RULE TO APPROXIMATE 

* THE MODIFIED MELLIN TRANSFORM. THE MODIFICATION 

* IS FREQUENCE TIMES THE FREQUENCY DERIVATIVE OF 

* THE FFT. 



******************************************************* 



SUBROUTINE SIMP (H, N ETS , NCOEF) 

DIMENSION H (NPTS) ,ID(5) 

COMMON XFM (256 ,2) ,CPHI (256, 128) ,SPHI (256, 128) , PI, 
1IXT (10) ,IYl (10) .KEY 

DATA ID/' I-flE ' , ' LLI N' , ' FRE • , 'QUEN* , • CY •/ 

CALL TITLE (ID) 



N 1 = NPTS -1 
N2 = NPTS - 2 



DO 100 J=1, NCOEF 
FSIMP = 2.0 

CMEGA = PI * FLOAT (J) / 18. 

COLOR = 1 . 

XFM (J, 1) = Q. 

XFM (J,2) = 0. 

DC 200 1*1,81 
I Ml = I 
10 = 1 + I 

IF (FSIMP .GT. 3.) GO- TC 67 
FSIMP = 4. 

GO TO 69 

67 FSIMP = 2. 

69 CONTINUE 

DELTA = H (10) - HfIMI) 

XFM ( J , 1 ) = XFM ( J # 1 ) + FSIMP * DELTA * CPHI(J,IM1) 

XFM (J, 2) = XFM(J,2) + FSIl^P * DELTA * S?HI(J,IM1) 

200 CONTINUE 

100 CONTINUE 
SETURN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 



** ******* ******* ************************************ 

* THE X A3 SUEBOUTINE TAKES THE MAGNITUDE OF THE 

* COMMON COMPLEX ABRAY XFM (256,2) ABC PLACES 

* THESE VALUES IN OUTPUT VECTOB <XHAG> FOB LATEfi 

* PRINTING. 



SUBROUTINE XAB(XMAG,NPTS) 

DIMENSION XMAG(NPTS) 

COMMON XFM (256. 2) ,CPHI (256,128) , SPHI (256, 128) ,PI, 
II XT (10) ,IYT (10f ,K&Y 
DO TOO I=1,NPTS 

XMAG(I) = 5QBT (XFM (1,1) **2+XFM (1 , 2J**2) 

100 CONTINUE 
BETUBN 
END 
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c 

C ******************************************************** 

C * THE STOW SUBROUTINE STORES THE INPUT ARRAY PRT(NPTS) 

C * INTO THE LOGICAL DEVICE 2 FOR LATER USE IN PLOTTING. 

C * NPLOT NUMEERS THE PLOTS FOR LATER IDENTIFICETION , 

C * AND A PLOT TITLE FOR THE MELLIN FREQUENCY IS ADDED 

C * FCR CONVENIENCE. 

C * THE MODULUS OF THE 

C * TRANSFORM TAKEN. SCALED TO UNIT MAGNITUDE. AND 
C * CUTPUT TO LOGICAL UNIT 2. 

C * INPUT: PRT - TO BE SCALED TO 1 AND WRITTEN 

C * TO LOGICAL UNIT 2. 

C * NPLOT - THE NUMBER OF THE PLOT 

C * KEY - NUMBER OF CURVE THIS PLOT 

C ******************************************************** 

c 



SUBROUTINE STOW (PRT ,NPTS) 

DIMENSION PRT (NPTS) 

COMMON XFM {256,2) , CPHI ( 256, 1 28) , SPHI (256,128) ,PI, 
1IIT (10) ,IYT (10) /KEY 





WRITE (3,13) 


KEY 


CKK 


IF (KEY . NE 


1) GO TO 12 


C 


WRITE AXIS 


LABELS 




WRITE (3,10) 
WRITE j3 ,10) 
FORMAT (10A4 


(IXT(I) ,1=1,10) 
(IYT(I) ,1=1,10) 


10 




12 


NPLOT = NP 
BPRT =0.0 


li0T + 1 


13 


FORMAT (14) 
DO 100 1=1, 


NPTS 


100 


IF (PRT (I) 
CONTINUE 


.GT. BPRT) BPRT 



20 

200 



IF (BPRT .LT. .00001) BPRT 
DO 200 1=1 , NPTS 
PBT(I) = PRT (I) / BPRT 
T = PLOAT(I) 

WRITE (3,20)1, PRT (I) 

FORMAT (2 (3X,F9.5) ) 

CONTINUE 

RETURN 

END 



PRT (I) 

1 . 
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c 

c 

c *********************************************** 

C * THE HOLD SUBROUTINE TAKES THE INPUT FILE * 

C * <FILIN> AND STORES IT IN THE OUTPUT FILE * 

C * <FILOUT> FOR TEMPORARY STORAGE. T BE FILE * 

C * <FILIN> REMAINS UNCHANGEE. * 

C ********************************** ******* ****** 

c 



SUBROUTINE HOLD (FILIN , FILOUT , NPTS) 
DIMENSION FILIN (NPTS) , FILOUT (NPTS) 



DO IOC I = 1, NPTS 
FILOUT (I) = FILIN (I) 

100 CONTINUE 
RETURN 
END 



C 

C 

C 

C 

C 

C 

C 

c 

c 



****************************************************** 



* THE ALTER SUBROUTINE MILL ALTER THE COMMON ARRAY 

* <WFM> AS SPECIFIED BY THE INPUT VARIABLES 

* <SCALE> AND <SHIFT> . THE ALTERED WFM IS OUTPUT 

* IN THE VECTOR <ALT> . 



****************************************************** 



100 



SUBROUTINE ALTER (ALI, MF M , SCALE, SHIFT, NPTS) 
DIMENSION ALT (NPTS) ,«FM (NPTS) .TOLD (256) .TNEii (256) 
COMMON XFM(256,2),CEHI (256, 128) ,SPHI (256,128) ,PI, 



1IXT (10) ,IYT (10 IMx 
DO 10O 1=1, NPTS 
IOLD(I) = FLOAT (I) 

TNEH(I) = TOLD (I) / SCALE + SHIFT 
CONTINUE 

CALL INTP3 (HFM, TOLD, ALT, TNEW, NPTS) 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 

c 

c 



** * ******* * * ****** * 7 *** * **** * * # **** * ** * A* * * * * * * * * * * ** * ** 

* INTP3 IS A SECOND OBDEE INTERPOLATION BASED OH A * 

* ECLINOMIAL EQUATING TO FIBST AND SECOND DERIVATIVES * 

* APPROXIMATED BY CENTRAL DIFFERENCES. THE INPUT * 

* VECTOR <XO> HAS OLD SAMPLES AT TIMES IN <TO>. THE * 

* NEW SAMPLE TIMES ARE INPUT THROUGH ARRAY <TN> AND * 

* THE COMPUTED SAMPLES AT THESE TIMES ARE OUTPUT IN * 

* ARE OUTPUT IN THE VECTOR <XN> AND <WFM>. * 

*** ************* *****$******** ****************$:********* 

SUBROUTINE INTP3 (XO.TO, XN,TN,NPTS) 

DIMENSION X3 (3) ,2o jf NPTS ), TO (NPTS) ,XN (NPTS) , TN (NPTS) 
COMMON XFM (256,2) ,CPHI (256, 128) ,5PHI (256, l28) , PI, 

1IXT (10) ,IYT (10) ,KEY 

CHOSE THE <IO> SAMPLE TIME CLOSES TO THE WARPED 
TIME HELD IN <TN> 

DC 60 1=1, NPTS 
DT = 100 



10 

25 



XPTS = FLOAT (NPTS) 



.AND. (TN (I) . LE. XPTS)) GO TO 5 



GO TO 60 
DO 10 J= 1 , NPTS 
DTABS = ABS (DT) 

CIST = TO (J). - TN (I) 
DIST = ABS (DIST) 

IF (DTABS .LT. DIST) 
DT = IN (I) - TO (J) 
CONTINUE 
JTIME = J - 1 
J1 = -1 
DO 30 J=1 ,3 
J 1 = J + JTIME- 1 



GO TO 25 





IF UJ1 . GE. 1) 
X3 (J) = 0.0 
GO TO 30 


29 


X3 (J) = XO (J1) 

CONTINUE 

TIMN = TN (I). 

1IMO = TO (JTIME 

XN (I) = FCN(X3, 

CONTINUE 

CO 500 1=1, NPTS 


30 


60 


500 


CONTINUE 

RETURN 

END 



.AND. (J1 .LE. NPTS) ) GO TO 29 



IN) 



12 S 



nnn non nnnn 



ECN IS THE INTERPOLAIIO N RULE 



FUNCTION FCN (X, TO , TN) 
DIMENSION X (3) 

COMPOTE THE COEFFICIENTS. 



A = (X (3) - 2.* X (2) + X (1 

/ 2.0 



A = (X (3) - 2.* X 



**2. - B 



2 



* A * 

* TO 



TO 



COMPUTE THE INTERPOLATED VALUE. 



FCN = A * (TN**2) + 3 * TN ♦ C 

BETUBN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



100 

200 

202 

101 



10 

20 



******************************************************** 

* THIS SUBROUTINE IS THE RESULT OF A CLOSED 

* FORM CALCULATION OF THE CONTINUOUS SINC**2 

* EEING TAKEN FROM THE TIME DOMAIN THROUGH THE 

* ENTIRE FOURIER- HELL IN PROCESSING. THE ALGORITHM 

* IS USED TO VERIFY THE FUM PROCESS. THE OUTPUT 

* FEATURE SPACE SHOULD BE IDENTICLE TO THAT PRODUCED 

* BY ANY OF THE MELLIN ALGORITHMS. FOR THIS REASON 

* THE SAMPLE POINTS ARE SYNCHRONIZED WITH THOSE USED 

* ABOVE. <NCOEF> FM COEFFICIENTS ARE USED. THE 

* COMMON VARIABLE <KEY> MUST BE SET TO 99 PRIOR TO 

* CALLING TO GET THE SINC**2 OUTPUT FEATURE SPACE. 

* TO OUTPUT THE CLOSED FORM SOLUTION FOR A RAMP IN 

* THE FREQUENCY DOMAIN. <KEY> MUST EQUAL 100. 
******************************************************** 



SUBROUTINE CFORM (CF , NCO EF) 

DIMENSION CE(NCOEF) 

COMMON XFM (256,2) , CPHI ( 256, 1 28) , SPHI (256,128) ,PI, 
1IXT (1C) ,IYT (10) , KEY 

IF (KEY .EQ. 100) GO TO 200 

IF (KEY .NE.99) RETURN 

DO 100 M= 1 , NCOEF 

XM = FLOAT (M) - 1.0 

HMAG = 1 + (PI * XM / 18.)**2 

CF(M) = SQRT(HMAG) / HMAG 

CONTINUE 

GO TO 101 

DO 202 I = 1, NCOEF 
XI = FLOAT (I) -1.0 
OMEGA = PI * XI / 18. 

CF (I) = 2 ./SQRT (4. + OMEGA**2) 

CONTINUE 

CONTINUE 

BEAD (2, 20) KEY 
READ (2 r 10 
READ (2, 10 
FORMAT (10 
FORMAT (14 

CALL STOW (CF, NCOEF) 

RETURN 
END 




127 



#«******«**** 



c 

c 

c 

c 

c 

c 



*** ***4** ********** **************V* ****** ******** 

* TITLE TITLES THE I AXIS FOR THE PLOTTING 

* EROGEAM ACCORDING TO THE ALGORITHM USED TO 

* GENERATE THE FEATURE SPACE. 

************************************************* 



SUBROUTINE TITLE (ID) 

DIMENSION ID (5) 

COMMON XFM (256,2) ,CPHI (256, 128) ,SPHI (256, 128) ,PI, 
1IXT MO) ,IYT j10) ,KEY " ' 

DO tOO 1*1-5 
IYT(I) = ID (I) 

100 CONTINUE 
RETURN 
END 
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